Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

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


From   Antoine Terracol <terracol@univ-paris1.fr>
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?

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/



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