Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.

# st: inaccurate hessian from mata: deriv()

 From Sun Yutao To Subject st: inaccurate hessian from mata: deriv() Date Thu, 25 Oct 2012 19:24:30 +0200

```Hello users,

I¡¯m having a problem on the function derive(), it looks like it¡¯s giving a
very inaccurate hessian. Here the problem goes:

I wrote a Stata probit log-like function and enveloped it with a Mata
function, which does no more than compute and sum the log-likelihood
contributions from the Stata program.
----------------------------------------------------------
function MleLikeMata(input_param,struct mPML_model scalar model,|LogLik) {
if (cols(input_param)!=1) {
param=input_param'
} else {
param=input_param
}

for(i=1;i<=rows(model.equation);i++) {
data_range=strtoreal(tokens(model.equation[i,7]))
param_range=model.location[i,1]::model.location[i,2]
if (model.equation[i,4]=="1") {
model.data[.,model.offset+i]=(model.data[.,data_range],J(rows(model.data),1,1))*param[param_range,.]
} else {			model.data[.,model.offset+i]=model.data[.,data_range]*param[param_range,.]

}
}
stata(model.StataLike+" "+model.XB)
LogLik=sum(model.data[.,model.offset])
//LogLik=model.data[.,model.offset]
if (args()==2) {
return(LogLik)
}
}
----------------------------------------------------------
Where model is just an extra argument I need

Then I use this code to compute the gradients and hessian:
----------------------------------------------------------
D=deriv_init()
deriv_init_evaluator(D,&MleLikeMata())
//deriv_init_evaluatortype(D,"v")
deriv_init_search(D, "off")deriv_init_params(D,
param')
deriv_init_argument(D, 1, model)
hess=deriv(D, 2)
----------------------------------------------------------

It looks like it's giving a more or less OK gradient but a VERY inaccurate
hessian (I'm comparing this with stata ml command with the same data and same
parameter vector, but -ml- works well!) , i.e. I'm 100% sure the Stata probit
log-likelihood function is concave and hence the hessian should be
negative-definite.
however, it gives some hessian like this.
+----------------------------------------------------------------------------+
1 |   238418.5803
2 |  -119209.2902    476837.1606
3 |  -178813.9352   -327825.5479    417232.5155
4 |  -59604.64508   -119209.2902   -208616.2578    357627.8705
5 |   149011.6127   -178813.9352   -59604.64508    29802.32254
0
+----------------------------------------------------------------------------+
The determinant is 2.43788e+27 which means it's certainly not
negative-definite (and you even have a 0 in the diagonal)

Possible fix I tried:

1. I tried to specify it as a v-type evaluator (one that returns a row vector)
2. I tried to scale sum of the log-likelihood by number of observations
3. change the 2nd argument in -deriv_init_search(D, "off")- to "bracket" or
"interpolate" (which gives even worse results)

Could anyone kindly explain me what's wrong and what I should do?

And infact I think I lost a lot features (fast convergence, etc.) from
Newton-Raphson due to the inaccurate hessian.

Oh in case you need, the Stata likelihood function is:
---------------------------------------------------
program define like, nclass
qui {
args lnf XB
replace `lnf' = ln( 1 - normprob(`XB'))     if \$ML_y1==0
replace `lnf' = ln(     normprob(`XB'))     if \$ML_y1==1
}
End
---------------------------------------------------
Which is the common one for probit

Best regards,
Sun Yutao

*
*   For searches and help try:
*   http://www.stata.com/help.cgi?search
*   http://www.stata.com/support/faqs/resources/statalist-faq/
*   http://www.ats.ucla.edu/stat/stata/
```