Laszlo Sandor <[email protected]> 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
[email protected]
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/
