[R] Hessian Matrix Issue

[dave fournier]

I wonder if your code is correct?

I ran your script until an error was reported. the data set
of 30 obs was


[1]  0  0  1  3  3  3  4  4  4  4  5  5  5  5  5  7  7  7  7  7  7  8  9 
10 11
[26] 12 12 12 15 16

I created a tiny AD Model Builder program to do MLE on it.

DATA_SECTION
   init_int nobs
   init_vector y(1,nobs)
PARAMETER_SECTION
   init_number log_mu
   init_number log_alpha
   sdreport_number mu
   sdreport_number tau
   objective_function_value f
PROCEDURE_SECTION
   mu=exp(log_mu);
   tau=1.0+exp(log_alpha);
   for (int i=1;i<=nobs;i++)
   {
     f-=log_negbinomial_density(y(i),mu,tau);
   }
It converged quickly and

The eigenvalues of the Hessian were

     4.711089774    78.27632341

and the estimates and std devs of the parameters mu and tau were

index   name       value      std dev

      3   mu         6.6000e+00 7.7318e-01
      4   tau        2.7173e+00 7.8944e-01

where tau is the variance divided by the mean.

This was all so simple that I suspect your (rather difficult to read)
R code is wrong, otherwise R must really suck at this kind of problem.

       Dave