st: [Mata newbie] optimize() vs -ml-

 From Antoine Terracol To statalist@hsphsun2.harvard.edu Subject st: [Mata newbie] optimize() vs -ml- Date Tue, 16 Oct 2007 20:17:28 +0200

```Dear _all,

I was wondering if there was any general rule to decide wether a given
likelihood maximisation problem will be more efficiently handled using
Mata's optimize() or Stata's -ml- command ?

When I compare the example code given in the help file for optimize() to
fit a ml estimation of a beta density (v0 method) whith the equivalent
-ml- code (lf method), I find that optimize() runs approximately 3 to 5
times faster than -ml-, depending on the size of the (simulated) dataset.

However, when fitting a simple probit model (two explanatory variables),
optimize() (v0 method) runs significantly slower than -ml- (lf method).

Is it caused by the matrix multiplication (cross()) needed to calculate
x*beta when explanatory variables are present, or is my code badly
written?

Antoine

---------Beta example------------------------------
. clear

. clear mata

. set obs 2000
obs was 0, now 2000

. g p=uniform()

. gen beta=invibeta(2,3,p)

. timer clear 1

. timer clear 2

. timer on 1

. mata:
------------------------------------------------- mata (type end to
exit) --------------------------------------------------
: mata clear

: void lnbetaden0(todo, beta,  data ,  lnf, S, H)
>  {
>          a   = beta[1]
>          b   = beta[2]
>          lnf = lngamma(a+b) :- lngamma(a) :- lngamma(b) :+
>                (a-1)*log(data) :+ (b-1)*log(1:-data)
>  }
note: argument todo unused
note: argument S unused
note: argument H unused

: x=0

: st_view(x,  ., "beta")

: S = optimize_init()

: optimize_init_evaluatortype(S, "v0")

: optimize_init_params(S, (0.1,0.5))

: optimize_init_argument(S, 1, x)

: betahat = optimize(S)
Iteration 0:  f(p) = -2320.5911
Iteration 1:  f(p) =  429.95828
Iteration 2:  f(p) =  470.22217
Iteration 3:  f(p) =  470.28914
Iteration 4:  f(p) =  470.28915

: sigmahat=sqrt(diagonal(optimize_result_V_oim(S)))'

: tstat=betahat:/sigmahat

: pval=2:*(1:-normal(abs(tstat)))

: result=betahat\sigmahat\tstat\pval

: result'
1             2             3             4
+---------------------------------------------------------+
1 |  2.001305167   .0594920662   33.63986652             0  |
2 |  3.002087416    .092931678   32.30424199             0  |
+---------------------------------------------------------+

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

. timer off 1

.
. cap prog drop mlbeta

. program define mlbeta
1. args lnf alpha beta
2. qui replace
lnf'=lngamma(`alpha'+`beta')-lngamma(`alpha')-lngamma(`beta')+(`alpha'-1)*log(\$ML_y1)+(`beta'-1)*log(1-\$ML_y1)
3. end

. ml model lf mlbeta (alpha: beta=) / beta

. mat I= (0.1 , 0.5)

. ml init I, copy

. timer on 2

. ml max, search(off) nolog

Number of obs   =
2000
Wald chi2(0)    =
.
Log likelihood =  470.28915                       Prob > chi2     =
.

------------------------------------------------------------------------------
|      Coef.   Std. Err.      z    P>|z|     [95% Conf.
Interval]
-------------+----------------------------------------------------------------
alpha        |
_cons |   2.001305   .0594921    33.64   0.000     1.884703
2.117907
-------------+----------------------------------------------------------------
beta         |
_cons |   3.002087   .0929317    32.30   0.000     2.819945
3.18423
------------------------------------------------------------------------------

. timer off 2

.
. timer list 1
1:      0.17 /        1 =       0.1720

. timer list 2
2:      0.64 /        1 =       0.6410

-----Probit example----------------------------------------

. clear

. clear mata

. set obs 2000
obs was 0, now 2000

. drawnorm x1 x2 eps

. g cons=1

. gen y=(1+x1+x2+eps>0)

. timer clear 1

. timer clear 2

. timer on 1

. mata:
------------------------------------------------- mata (type end to
exit) --------------------------------------------------
: void lnprob0(todo, beta, y, x , lnf, S, H)
>  {
>                 xb=cross(x', beta')
>                 lnf = ln(normal(xb)):*y:+ln(normal(-xb)):*(1:-y)
>  }
note: argument todo unused
note: argument S unused
note: argument H unused

: y=x=.

: st_view(y,  ., "y")

: st_view(x,  ., ("x1","x2","cons"))

: S = optimize_init()

: optimize_init_evaluator(S, &lnprob0())

: optimize_init_evaluatortype(S, "v0")

: optimize_init_params(S, (0.1,0.5,0))

: optimize_init_argument(S, 1, y)

: optimize_init_argument(S, 2, x)

: betahat = optimize(S)
Iteration 0:  f(p) = -1149.5343
Iteration 1:  f(p) = -741.17206
Iteration 2:  f(p) = -723.61705
Iteration 3:  f(p) = -723.55496
Iteration 4:  f(p) = -723.55495

: sigmahat=sqrt(diagonal(optimize_result_V_oim(S)))'

: tstat=betahat:/sigmahat

: pval=2:*(1:-normal(abs(tstat)))

: (betahat\sigmahat\tstat\pval)'
1             2             3             4
+---------------------------------------------------------+
1 |  .9336980122   .0492557892   18.95610703             0  |
2 |  .9553768942   .0485996012   19.65812208             0  |
3 |   .933645069   .0443509123   21.05131597             0  |
+---------------------------------------------------------+

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

. timer off 1

. capture program drop myprobit

. program define myprobit
1. args lnf xb
2. quietly replace `lnf'=ln(normprob(`xb'))*y+ln(normprob(-`xb'))*(1-y)
3. end

. ml model lf myprobit (y=x1 x2)

. mat I= (0.1 , 0.5,0)

. ml init I, copy

. timer on 2

. ml max, search(off) nolog

Number of obs   =
2000
Wald chi2(2)    =
534.52
Log likelihood = -723.55495                       Prob > chi2     =
0.0000

------------------------------------------------------------------------------
y |      Coef.   Std. Err.      z    P>|z|     [95% Conf.
Interval]
-------------+----------------------------------------------------------------
x1 |    .933698   .0492558    18.96   0.000     .8371584
1.030238
x2 |   .9553769   .0485996    19.66   0.000     .8601234
1.05063
_cons |   .9336451   .0443509    21.05   0.000     .8467189
1.020571
------------------------------------------------------------------------------

. timer off 2

.
.
. timer list 1
1:      0.69 /        1 =       0.6880

. timer list 2
2:      0.36 /        1 =       0.3600

*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/
```