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?
Thanks in advance,
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_evaluator(S, &lnbetaden0())
: 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/
