Re: st: question concerning translation from Matlab to Mata

[DC <dcase79@gmail.com>]

Hi Matthew,

I'll check this out. Thanks a lot for the help.

Regards,
Marcus

On Fri, Mar 11, 2011 at 9:43 AM, Matthew J Baker
<matthew.baker@hunter.cuny.edu> wrote:
> After trouble shooting the program, I'm not sure exactly what I
> did that fixed it, but I got the following example to run. One
> additional problem you were having is that the initial parameter
> vector has to be 2*cols(X) long rather than just cols(X) long.
> More conceptually, it looks like you wrote
> d=1+exp(xb1)+exp(xb1) instead of the correct expression
> d=1+exp(xb1)+exp(xb2). But that won't stop the program from
> running. In any case, here's the example I got to work (using
> Cameron and Trivedi's data), which reproduces the results of
> the mlogit command in stata:
>
> clear all
> use "http://cameron.econ.ucdavis.edu/musbook/mus15data.dta";
> tab mode, gen(dum)
> gen y=0
> replace y=1 if dum1==1
> replace y=2 if dum2==1
> mlogit y income
> mata
> st_view(y=.,.,"y")
> st_view(x=.,.,"income")
> x=x,J(rows(y),1,1)
> void lmle(todo,b,y,x,lnf,S,H)
> {
>        k=cols(x)
>        b1=b[1,1::k]
>        b2=b[1,k+1::2*k]
>        xb1=x*b1'
>        xb2=x*b2'
>        d=1:+exp(xb1)+exp(xb2)
>        p_0 = 1:/d
>        p_1=exp(xb1):/d
>        p_2=exp(xb2):/d
>        lnf=(y:==0):* log(p_0):+(y:==1):*log(p_1):+
> (y:==2):*log(p_2)
> }
> S=optimize_init()
> optimize_init_evaluator(S,&lmle())
> optimize_init_evaluatortype(S,"gf0")
> optimize_init_argument(S,1,y)
> optimize_init_argument(S,2,x)
> optimize_init_params(S,J(1,2*cols(x),0))
> b=optimize(S)
> st_matrix("b",b)
> end
> ereturn post b
> ereturn display
>
> Hope that helps!
>
> Matt Baker
>
> Dr. Matthew J. Baker
> Department of Economics
> Hunter College and the Graduate Center, CUNY
>
> ---- Original message ----
>>Date: Fri, 11 Mar 2011 08:51:16 -0500
>>From: owner-statalist@hsphsun2.harvard.edu (on behalf of DC
> <dcase79@gmail.com>)
>>Subject: Re: st: question concerning translation from Matlab to
> Mata
>>To: statalist@hsphsun2.harvard.edu
>>
>>Hi Brian and Matthew,
>>
>>Thanks a lot for the help. I tried those suggestions and
> unfortunately
>>I am still getting the following error:
>>
>> logitmle():  3301  subscript invalid
>>
>>Thanks,
>>Marcus
>>
>>On Thu, Mar 10, 2011 at 10:39 PM, Matthew J Baker
>><matthew.baker@hunter.cuny.edu> wrote:
>>> Another problem that might produce the conformability error
> is the use of (y==0). I'm guessing that this might be better
> coded as (y:==0).
>>>
>>> MJB
>>>
>>> Dr. Matthew J. Baker
>>> Department of Economics
>>> Hunter College and the Graduate Center, CUNY
>>>
>>>
>>> ---- Original message ----
>>>>Date: Thu, 10 Mar 2011 18:09:53 -0500
>>>>From: owner-statalist@hsphsun2.harvard.edu (on behalf of
> "Brian P. Poi" <brian@poiholdings.com>)
>>>>Subject: Re: st: question concerning translation from Matlab
> to Mata
>>>>To: statalist@hsphsun2.harvard.edu
>>>>
>>>>
>>>>
>>>>On 3/10/2011 5:55 PM, DC wrote:
>>>>> Hi All,
>>>>>
>>>>> I'm attempting to translate some code for an optimal
> stopping problem
>>>>> that I had previously written in Matlab to Mata and I've
> encountered a
>>>>> problem
>>>>> in writing down a likelihood function.
>>>>>
>>>>> A simplified version in Matlab would be a multinomial logit
> with three
>>>>> choices, choice 0 parameter normalized to zero, :
>>>>> ...
>>>>> [n k]=size(x);
>>>>> D = 1+exp(x*b(1:k))+exp(x*b(k+1:2*k));
>>>>> p_0 = 1./D;
>>>>> p_1 = exp(x*b(1:k))./D;
>>>>> p_2 = exp(x*b(k+1:2*k))./D;
>>>>> lnlike=-sum( (y==0).*log(p_0) + (y==1).*log(p_1) +
> (y==2).*log(p_2));
>>>>> ....
>>>>> However is it not clear how to specify and estimate
> something like
>>>>> this using Mata
>>>>> I tried the following:
>>>>>
>>>>> void lmle(todo,b,y, x,lnf, S, H)
>>>>> {
>>>>>      k = cols(x)
>>>>>      b1 = b[1::k,1]
>>>>>      b2 = b[k+1::2*k,1]
>>>>>      xb1 = x*b1'
>>>>>      xb2 = x*b2'
>>>>>
>>>>>      d = 1 :+ (exp(xb1) + exp(xb1))
>>>>>      p_0 = 1 :/ d
>>>>>      p_1 = exp(xb1) :/ d
>>>>>      p_2 = exp(xb2) :/ d
>>>>>      lnf = (y==0) :* log(p_0) + (y==1) :* log(p_1)  +
> (y==2) :* log(p_2)
>>>>> }
>>>>>
>>>>> S = optimize_init()
>>>>> optimize_init_evaluator(S,&lmle())
>>>>> optimize_init_evaluatortype(S, "gf0")
>>>>> optimize_init_argument(S,1,y)
>>>>> optimize_init_argument(S,2,x)
>>>>> optimize_init_params(S, J(1,cols(x),0))
>>>>>
>>>>> b = optimize(S)
>>>>>
>>>>
>>>>
>>>>One immediate problem is that -optimize- passes the
> parameters as a row
>>>>vector, not a column vector, so the lines
>>>>
>>>> >     b1 = b[1::k,1]
>>>> >     b2 = b[k+1::2*k,1]
>>>>
>>>>should probably be
>>>>
>>>>       b1 = b[1, 1::k]
>>>>       b2 = b[1, k+1::2*k]
>>>>
>>>>That may or may not be the only problem.
>>>>
>>>>    -- Brian Poi
>>>>    -- brian@poiholdings.com
>>>>
>>>>*
>>>>*   For searches and help try:
>>>>*   http://www.stata.com/help.cgi?search
>>>>*   http://www.stata.com/support/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/statalist/faq
>>> *   http://www.ats.ucla.edu/stat/stata/
>>>
>>
>>
>>
>>--
>>Marcus Casey, Ph.D.
>>Duke University
>>
>>*
>>*   For searches and help try:
>>*   http://www.stata.com/help.cgi?search
>>*   http://www.stata.com/support/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/statalist/faq
> *   http://www.ats.ucla.edu/stat/stata/
>



-- 
Marcus Casey, Ph.D.
Duke University

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