st: inaccurate hessian from mata: deriv()

[Sun Yutao <yutao.sun.statalist@outlook.com>]

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)
grad=deriv(D, 1)'
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/