Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: st: moptimize within moptimize?


From   jpitblado@stata.com (Jeff Pitblado, StataCorp LP)
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: moptimize within moptimize?
Date   Fri, 06 Nov 2009 14:32:22 -0600

Laszlo Sandor <sandorl@gmail.com> asks about using -moptimize()- within
an -moptimize()- evaluator:

> I am trying to code up a simple EM algorithm to estimate a latent
> class model (a.k.a finite mixture on a higher level than the
> observations). This should be doable if I can keep track of all my
> parameters in the score and Hessian (of simple classical normals) in
> the M-step (maximization of log-likelihood for a mixing taken as
> given). I'll try to keep the M-step in -moptimize-. But it would be
> also nice to embed the M-step within another moptimize call for the
> whole procedure. Do you think it is feasible, or it is unnecessary, or
> even unfeasible? It would be great to have some expert opinion on this
> before I go down that path.

Yes, you can nest calls to -moptimize()- within other -moptimize()-
evaluators.

Here is an example, where I use -moptimize()- to fit a linear regression,
treating the sigma parameter of the normal distribution as a nuisance
parameter by -moptimize()-ing the normal likelihood with respect to ln(sigma)
conditionally on the regression coefficients.

The Stata output from this example is included after my signature.

***** BEGIN:
sysuse auto

mata:

void eval1(M1, todo, b, v, g, H)
{
	y	= moptimize_util_depvar(M1)
	xb	= moptimize_util_xb(M1, b, 1)

	M2	= moptimize_init()
	moptimize_init_evaluator(M2, &eval2())
	moptimize_init_evaluatortype(M2, "d0")
	moptimize_init_depvar(M2, 1, y)
	moptimize_init_userinfo(M2, 1, xb)
	moptimize_init_eq_cons(M2, 1, "on")
	moptimize_init_tracelevel(M2, "none")

	if (_moptimize(M2)) {
		v = .
	}
	v	= moptimize_result_value(M2)
}

void eval2(M2, todo, b, v, g, H)
{
	y	= moptimize_init_depvar(M2, 1)
	xb	= moptimize_init_userinfo(M2, 1)
	lnsigma	= moptimize_util_xb(M2, b, 1)

	v	= moptimize_util_sum(M2, lnnormalden(y, xb, exp(lnsigma)))
}

end

ml model d0 eval1() (mpg = turn trunk displ), maximize
ml display

regress mpg turn trunk displ

***** END:

--Jeff
jpitblado@stata.com

Stata output from the preceeding example:

. sysuse auto
(1978 Automobile Data)

. 
. mata:
------------------------------------------------- mata (type end to exit) -----
: 
: void eval1(M1, todo, b, v, g, H)
> {
>         y       = moptimize_util_depvar(M1)
>         xb      = moptimize_util_xb(M1, b, 1)
> 
>         M2      = moptimize_init()
>         moptimize_init_evaluator(M2, &eval2())
>         moptimize_init_evaluatortype(M2, "d0")
>         moptimize_init_depvar(M2, 1, y)
>         moptimize_init_userinfo(M2, 1, xb)
>         moptimize_init_eq_cons(M2, 1, "on")
>         moptimize_init_tracelevel(M2, "none")
> 
>         if (_moptimize(M2)) {
>                 v = .
>         }
>         v       = moptimize_result_value(M2)
> }
note: argument todo unused
note: argument g unused
note: argument H unused

: 
: void eval2(M2, todo, b, v, g, H)
> {
>         y       = moptimize_init_depvar(M2, 1)
>         xb      = moptimize_init_userinfo(M2, 1)
>         lnsigma = moptimize_util_xb(M2, b, 1)
> 
>         v       = moptimize_util_sum(M2, lnnormalden(y, xb, exp(lnsigma)))
> }
note: argument todo unused
note: argument g unused
note: argument H unused

: 
: end
-------------------------------------------------------------------------------

. 
. ml model d0 eval1() (mpg = turn trunk displ), maximize

initial:       log likelihood = -333.93641
alternative:   log likelihood = -332.30036
rescale:       log likelihood = -257.15293
Iteration 0:   log likelihood = -257.15293  (not concave)
Iteration 1:   log likelihood = -227.71457  (not concave)
Iteration 2:   log likelihood = -221.58427  
Iteration 3:   log likelihood = -210.44683  
Iteration 4:   log likelihood = -207.18345  
Iteration 5:   log likelihood = -202.23753  
Iteration 6:   log likelihood = -201.61888  
Iteration 7:   log likelihood = -201.61776  
Iteration 8:   log likelihood = -201.61776  

. ml display

                                                  Number of obs   =         74
                                                  Wald chi2(3)    =     105.46
Log likelihood = -201.61776                       Prob > chi2     =     0.0000

------------------------------------------------------------------------------
         mpg |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        turn |  -.4971904   .1612735    -3.08   0.002    -.8132806   -.1811003
       trunk |   -.222583   .1316643    -1.69   0.091    -.4806403    .0354744
displacement |  -.0196434   .0077819    -2.52   0.012    -.0348956   -.0043913
       _cons |   47.94784   5.145066     9.32   0.000      37.8637    58.03199
------------------------------------------------------------------------------

. 
. regress mpg turn trunk displ

      Source |       SS       df       MS              Number of obs =      74
-------------+------------------------------           F(  3,    70) =   33.25
       Model |  1435.86938     3  478.623127           Prob > F      =  0.0000
    Residual |  1007.59008    70   14.394144           R-squared     =  0.5876
-------------+------------------------------           Adj R-squared =  0.5700
       Total |  2443.45946    73  33.4720474           Root MSE      =   3.794

------------------------------------------------------------------------------
         mpg |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        turn |  -.4971904    .165826    -3.00   0.004      -.82792   -.1664608
       trunk |   -.222583   .1353748    -1.64   0.105    -.4925796    .0474136
displacement |  -.0196434   .0080013    -2.46   0.017    -.0356016   -.0036853
       _cons |   47.94784   5.290301     9.06   0.000     37.39667    58.49902
------------------------------------------------------------------------------

<end>
*
*   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/



© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index