Re: st: AW: Logit maximum likelihood in mata code

[Aljar Meesters <aljar.meesters@gmail.com>]

Dear Nick,

There are two issues that bite you. One is that you have to extend the
initial parameter vector with one extra element, since you are
estimating two parameters. If you fix this you will get a
conformability error, since you treat b as a column vector while it is
a row vector.
This should work
----
clear all
sysuse auto, clear
gen one = 1
mata
st_view(X=0,.,("weight", "one"))
st_view(y=0,.,("foreign"))
void logistic(todo, p, y, X, lf, g, H)
{
b = p[1, (1::cols(X))]'    // Transpose such that b is a column vector
lf = y' * ln(invlogit(X*b)):+(1:-y)'*ln(1:-invlogit(X*b))
}
S = optimize_init()
optimize_init_evaluator(S, &logistic())
optimize_init_evaluatortype(S, "gf0")
optimize_init_params(S, (0, 0))     // One extra element. Starting
values are (0, 0)
optimize_init_argument(S, 1, y)
optimize_init_argument(S, 2, X)
optimize_init_which(S,"max")
p = optimize(S)
end
---





2013/4/19 nick bungy <nickbungystata@hotmail.co.uk>:
> Hi Klaus,
> Thanks for the reply. I intend to have a look through the support materials you linked. My aim here was to derive some basic MLEs in the optimizer to help get to grips with what it can do (as it seems a
> powerful tool)
> I cleaned up my code and the optimizer runs with just one indep variable (see first set of code) but the moment I try to add a second indep variable I just get invalid subscript. I am presumably making an error somewhere but struggling to find it.
> One variable - does work:
> sysuse auto, clear
> mata
> st_view(X=0,.,("weight"))
> st_view(y=0,.,("foreign"))
> void logistic(todo, p, y, X, lf, g, H)
> {
> b = p[1,1::1]
> lf = y'*ln(invlogit(X*b)):+(1:-y)'*ln(1:-invlogit(X*b))
> }
> S = optimize_init()
> optimize_init_evaluator(S, &logistic())
> optimize_init_evaluatortype(S, "gf0")
> optimize_init_params(S, 0)
> optimize_init_argument(S, 1, y)
> optimize_init_argument(S, 2, X)
> optimize_init_which(S,"max")
> p = optimize(S)
> --------------
> Two variables - subscript invalid:
> sysuse auto, clear
> gen one = 1
> mata
> st_view(X=0,.,("weight", "one"))
> st_view(y=0,.,("foreign"))
> void logistic(todo, p, y, X, lf, g, H)
> {
> b = p[1, (1::cols(X))]
> lf = y'*ln(invlogit(X*b)):+(1:-y)'*ln(1:-invlogit(X*b))
> }
> S = optimize_init()
> optimize_init_evaluator(S, &logistic())
> optimize_init_evaluatortype(S, "gf0")
> optimize_init_params(S, (0))
> optimize_init_argument(S, 1, y)
> optimize_init_argument(S, 2, X)
> optimize_init_which(S,"max")
> p = optimize(S)
>> From: kpforr@googlemail.com
>> To: statalist@hsphsun2.harvard.edu
>> Subject: st: AW: Logit maximum likelihood in mata code
>> Date: Fri, 19 Apr 2013 12:45:07 +0200
>>
>> <>
>> Dear Nicholas,
>>
>> I do not understand your code, but your problem is extensively discussed in
>> this Stata-Press-book:
>> http://www.statapress.com/books/maximum-likelihood-estimation-stata/
>> You find a working mata-code for your problem here.
>> http://www.statapress.com/data/ml4.html
>>
>> best
>>
>> Klaus
>>
>>
>> __________________________________
>>
>> Dr. Klaus Pforr
>> GESIS -- Leibniz Institut für Sozialwissenschaft
>> B2,1
>> Postfach 122155
>> D - 68072 Mannheim
>> Tel: +49 621 1246 298
>> Fax: +49 621 1246 100
>> E-Mail: klaus.pforr@gesis.org
>> __________________________________
>>
>> -----Ursprüngliche Nachricht-----
>> Von: owner-statalist@hsphsun2.harvard.edu
>> [mailto:owner-statalist@hsphsun2.harvard.edu] Im Auftrag von nick bungy
>> Gesendet: Freitag, 19. April 2013 09:34
>> An: statalist@hsphsun2.harvard.edu
>> Betreff: st: Logit maximum likelihood in mata code
>>
>> Dear Statalist,
>> I'm having trouble setting up the optimizer to find the maximum likelihood
>> of a logit model. I'm still learning the syntax and set up of mata, though
>> I've tried to follow 'help mf_optimize' as best I can.
>> Could someone assist in where I am going wrong with the code?
>> ----------
>> sysuse auto, clear
>> mata
>> st_view(x=0,.,("weight"))
>> st_view(y=0,.,("foreign"))
>> x=x,J(rows(x),1,1)
>> void logistic(todo, p, y, x, llf, g, H)
>> {
>> b = p[1,1::2]
>> One =J(rows(x),1,1)
>> llf = y'*log(One:/One+exp(-(x*b)))+(One-y)'*log(One:/One+exp(-(x*b)))
>> }
>> S = optimize_init()
>> optimize_init_evaluator(S, &logistic()) optimize_init_params(S, 0)
>> optimize_init_argument(S, 1, y) optimize_init_argument(S, 2, x)
>> optimize_init_which(S,"max")
>> p = optimize(S)
>> -----------
>> Many thanks,
>> Nicholas
>> *
>> * For searches and help try:
>> * http://www.stata.com/help.cgi?search
>> * http://www.stata.com/support/faqs/resources/statalist-faq/
>> * http://www.ats.ucla.edu/stat/stata/
>>
>>
>> *
>> * For searches and help try:
>> * http://www.stata.com/help.cgi?search
>> * http://www.stata.com/support/faqs/resources/statalist-faq/
>> * http://www.ats.ucla.edu/stat/stata/
> *
> *   For searches and help try:
> *   http://www.stata.com/help.cgi?search
> *   http://www.stata.com/support/faqs/resources/statalist-faq/
> *   http://www.ats.ucla.edu/stat/stata/

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/