Statalist


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

st: problem with MLE


From   DEEPANKAR BASU <basu.15@osu.edu>
To   statalist@hsphsun2.harvard.edu
Subject   st: problem with MLE
Date   Sun, 02 Dec 2007 17:20:00 -0500

I am trying to estimate a fertility model with birth history data using maximum likelihood. The data comes from a survey (DHS) and I am using Stata 8.1.

When I run *ml check*, I get the message that all tests have been passed. Then, once the optimization starts, the program spits out the error message that "numerical derivatives could not be calculated" and stops. Any suggestions on how to go about finding the error, correcting it and making the program work?

In Stata language, my likelihood evaluator is of the type *lf*; and so I have not programmed the derivative or the hessian. Below I give the likelihood evaluator program and also the log file containing the output.

Any suggestions to get this working will be most welcome.

Deepankar

--------- Start: Likelihood function -------------

clear
set more 1


capture program drop step1mle_lf
program step1mle_lf
   version 8.1
     

   #delimit ;
   args lnf theta1 theta2 theta3 theta4 theta5 theta6 theta7 theta8 theta9 
            theta10 theta11 theta12 theta13 theta14;
   
 /* theta4 = generates probability of male birth (will be estimated as a parameter) */

   
   tempvar q1 q2 q3 q4 q5 q6 q7 q8 q9 q10 q11 q12 q13 q14 x1 x2 x3 x4 x5 x6 
           p21 p31 p32 p41 p42 p43 p51 p52 p53 p54 p61 p62 p63 p64 p65;
   #delimit cr
     


   quietly {    /*--------------- BEGINNING OF QUIET COMPUTATION ---------------*/



/* ------------ GENERATING THE CONDITIONAL PROBABILITIES ------------------- */

      gen double `q1' = 1/(1+exp(`theta5'))   
      gen double `q2' = exp(`theta5')/(1+exp(`theta5'))   

      gen double `q3' = 1/(1+exp(`theta6')+exp(`theta7'))   
      gen double `q4' = exp(`theta6')/(1+exp(`theta6')+exp(`theta7'))   
      gen double `q5' = exp(`theta7')/(1+exp(`theta6')+exp(`theta7'))   

      gen double `q6' = 1/(1+exp(`theta8')+exp(`theta9')+exp(`theta10'))   
      gen double `q7' = exp(`theta8')/(1+exp(`theta8')+exp(`theta9')+exp(`theta10'))   
      gen double `q8' = exp(`theta9')/(1+exp(`theta8')+exp(`theta9')+exp(`theta10'))   
      gen double `q9' = exp(`theta10')/(1+exp(`theta8')+exp(`theta9')+exp(`theta10'))
      
      gen double `q10' = 1/(1+exp(`theta11')+exp(`theta12')+exp(`theta13')+ exp(`theta14'))   
      gen double `q11' = exp(`theta11')/(1+exp(`theta11')+exp(`theta12')+exp(`theta13')+ exp(`theta14'))   
      gen double `q12' = exp(`theta12')/(1+exp(`theta11')+exp(`theta12')+exp(`theta13')+ exp(`theta14'))   
      gen double `q13' = exp(`theta13')/(1+exp(`theta11')+exp(`theta12')+exp(`theta13')+ exp(`theta14'))   
      gen double `q14' = exp(`theta14')/(1+exp(`theta11')+exp(`theta12')+exp(`theta13')+ exp(`theta14'))   



/* -------------------------------------------------------------------------------- */
/* ----------------------GENERATING THE PROBABILITY TABLE  ------------------------ */
/* p`i'`k' IS THE CONDITIONAL PROBABILITY OF OBSERVING THE BIRTH SEQUENCE           */
/* GIVEN THAT THE DESIRED FAMILY SIZE IS `i' AND THE DESIRED TARGET FOR BOYS IS `k' */

/* STEP 1  */
   forval i=2/6 {
     local m = `i'-1
     forval j=1/`m' {
       gen double `p`i'`j'' = 0
     }
   }


/* STEP 2   */
   tempvar cboys c1 c2 c3 c4 c5 c6
   gen byte `cboys' = 0

   gen `c1' = ($ML_y5==1)
   replace `cboys' = `cboys' + `c1'
   replace `c1' = `c1'*`cboys'

   gen `c2' = ($ML_y6==1)
   replace `cboys' = `cboys' + `c2'
   replace `c2' = `c2'*`cboys'

   gen `c3' = ($ML_y7==1)
   replace `cboys' = `cboys' + `c3'
   replace `c3' = `c3'*`cboys'

   gen `c4' = ($ML_y8==1)
   replace `cboys' = `cboys' + `c4'
   replace `c4' = `c4'*`cboys'

   gen `c5' = ($ML_y9==1)
   replace `cboys' = `cboys' + `c5'
   replace `c5' = `c5'*`cboys'

   gen `c6' = ($ML_y10==1)
   replace `cboys' = `cboys' + `c6'
   replace `c6' = `c6'*`cboys'


/*---- PROBABILITY OF BIRTH OF A MALE OFFSPRING ----*/
   tempvar pmale
   gen `pmale' = 1/(1+exp(`theta4'))


   forval i=1/6 {
       forval j=1/6 {
        forval l=1/6 {  
         local m=`i'-1
          forval k=1/`m' {
             replace `p`i'`k''=((`pmale')^`j')*((1-`pmale')^`l') if ($ML_y1==`i' & $ML_y4==`j' & $ML_y11==`l')
          }
        }
     }
   }


/* STEP 3  */
   forval i=1/6 {
     forval j=1/6 {
       forval l=1/6 {
          local m = `i'-1
          forval k=1/`m' {
             replace `p`i'`k''=0             if ($ML_y1==`i' & $ML_y2==`j' & $ML_y1==$ML_y2 & `c`j''!=`k' & $ML_y4>=`k')
             replace `p`i'`k''=0             if ($ML_y1==`i' & $ML_y2==`j' & $ML_y1>$ML_y2 & `c`j''!=`k')
             replace `p`i'`k''=((`pmale')^`j')*((1-`pmale')^`l') if ($ML_y1==`i' & $ML_y4==`j' & $ML_y11==`l' & 

$ML_y1==$ML_y2 & $ML_y2==0)
          }
        }
     }
   }

/* ----------------- END OF THE PART THAT GENERATES THE PROBABILITY TABLE --------------------------- */



      gen double `x1' = (`pmale'^($ML_y4))*((1-`pmale')^($ML_y11))
      gen double `x2' = exp(lnfact($ML_y1))
      gen double `x3' = `x1'*`x2'     
      gen double `x4' = ((exp(`theta1'))^($ML_y1))/exp(exp(`theta1'))
      gen double `x5' = norm(-`theta3')
      gen double `x6' = ((exp(`theta1'+`theta2'))^($ML_y1))/exp(exp(`theta1'+`theta2'))
      


      }  /* ----------------- END OF QUIET COMPUTATION ---------------- */


      #delimit ;
      quietly replace `lnf' = ln( ($ML_y3*(`x1')*(1/`x2')*`x4'*`x5') + ((1-`x5')*`x6'*(1/`x2')*( `p21' + `p31'*`q1'           

                              + `p32'*`q2' + `p41'*`q3' + `p42'*`q4' + `p43'*`q5' + `p51'*`q6' + `p52'*`q7'
                              + `p53'*`q8' + `p54'*`q9' + `p61'*`q10' + `p62'*`q11' + `p63'*`q12'
                              + `p64'*`q13' + `p65'*`q14')));
      #delimit cr
   
end  

---------- End: Likelihood function ------------- 


---------- Start: Log file ----------------------
       log:  c:\deepankar\gender-data\step1\india-mle-step1-results-1998.smcl
  log type:  smcl
 opened on:   2 Dec 2007, 16:55:03

. 
. set mem 100m

Current memory allocation

                    current                                 memory usage
    settable          value     description                 (1M = 1024k)
    --------------------------------------------------------------------
    set maxvar         5000     max. variables allowed           1.733M
    set memory          100M    max. data space                100.000M
    set matsize         400     max. RHS vars in models          1.254M
                                                            -----------
                                                               102.987M

. use india-step1-mle-data-1998.dta   /* DATSET: Contains observations with N=2,3,4,5,6 */

. 
. 
. /*---------- DETAILS OF SURVEY DESIGN -------------*/
. gen smpwt = v005/1000000

. svyset [pweight=smpwt]
pweight is smpwt

. 
. 
. /* ------------------ UNRESTRICTED MODEL ---------------------------*/
. #delimit ;
delimiter now ;
. ml model lf step1mle_lf (fsize: dfsize = edu age edu rur work middle poor hindu, nocons) (alpha:) 
> (target: alive dfsize1 nboy cord1 cord2 cord3 cord4 cord5 cord6 ngirl = age edu rur work middle poor contra jfam hindu
> ) 
>  /four /five /six /seven /eight /nine /ten /eleven /twelve /thirteen /fourteen, svy robust tech(bhhh nr);
note: edu dropped due to collinearity

.   #delimit cr
delimiter now cr
. 
. ml check

Test 1:  Calling step1mle_lf to check if it computes log pseudolikelihood and
         does not alter coefficient vector...
         Passed.

Test 2:  Calling step1mle_lf again to check if the same log pseudolikelihood value
         is returned...
         Passed.

Test 3:  Calling step1mle_lf to check if 1st derivatives are computed...
         test not relevant for method lf.

Test 4:  Calling step1mle_lf again to check if the same 1st derivatives are
         returned...
         test not relevant for method lf.

Test 5:  Calling step1mle_lf to check if 2nd derivatives are computed...
         test not relevant for method lf.

Test 6:  Calling step1mle_lf again to check if the same 2nd derivatives are
         returned...
         test not relevant for method lf.

------------------------------------------------------------------------------
Searching for alternate values for the coefficient vector to verify that
step1mle_lf returns different results when fed a different coefficient vector:

Searching...
initial:       log pseudolikelihood =     -<inf>  (could not be evaluated)
searching for feasible values +

feasible:      log pseudolikelihood = -29212.563
improving initial values ..........
improve:       log pseudolikelihood = -29212.563

continuing with tests...
------------------------------------------------------------------------------

Test 7:  Calling step1mle_lf to check log pseudolikelihood at the new values...
         Passed.

Test 8:  Calling step1mle_lf requesting 1st derivatives at the new values...
         test not relevant for method lf.

Test 9:  Calling step1mle_lf requesting 2nd derivatives at the new values...
         test not relevant for method lf.

------------------------------------------------------------------------------
                         step1mle_lf HAS PASSED ALL TESTS
------------------------------------------------------------------------------

Test 10: Does step1mle_lf produce unanticipated output?
         This is a minor issue.  Stata has been running step1mle_lf with all
         output suppressed.  This time Stata will not suppress the output.
         If you see any unanticipated output, you need to place quietly in
         front of some of the commands in step1mle_lf.

-------------------------------------------------------------- begin execution
---------------------------------------------------------------- end execution

. 
. ml search
initial:       log pseudolikelihood = -29212.563
improve:       log pseudolikelihood = -26594.685
rescale:       log pseudolikelihood = -17651.076
rescale eq:    log pseudolikelihood = -15128.187

. 
. ml max, difficult

initial:       log pseudolikelihood = -15128.187
rescale:       log pseudolikelihood = -15128.187
rescale eq:    log pseudolikelihood = -15128.187
(setting optimization to BHHH)
could not calculate numerical derivatives
missing values encountered
r(430);

end of do-file
r(430);

. log close
       log:  c:\deepankar\gender-data\step1\india-mle-step1-results-1998.smcl
  log type:  smcl
 closed on:   2 Dec 2007, 17:04:00

-------------- End: Log file ----------------



*
*   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–2021 StataCorp LLC   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index