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
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.