Dear Filipe On 14 October 2013 17:42, Filipe Ribeiro <flipjribe...@hotmail.com> wrote: > Dear Arne, > > First of all, thank you very much for your reply. Secondly, please accept my > apologies for only give you an answer now, but I was moving from one country > to another and only today I was able to get back to work. > I did not write the model specification because it was just a trial, but my > main idea was to apply the "Nelder-Mead" model. > Since the beginning my guess was that something is wrong with the > log-likelihood function, however, I don't find any problem with the same > exact function applying the "optim" function: > > loglike.GGompi <- function(theta,age,deaths,exposures) { > alpha <- exp(theta[1]) > beta <- exp(theta[2]) > gamma <- exp(theta[3]) > first <- alpha*exp(beta*age) > second <- 1+(((alpha*gamma)/beta)*(exp(beta*age)-1)) > mu <- first/second > llk <- -sum((deaths * log(mu)) + (- mu*exposures)) > return(llk) > } > > > fit1 <- optim(par=c(-4.1402817, -0.6375773, -1.6945914), > loglike.GGompi, > age=0:6, > deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198), > exposures=c(935023.67, 819392.00, 724568.17, 470947.00, > 231951.64, 69502.65, 15798.72)) > > exp(fit1$par) > >> exp(fit1$par) > [1] 0.01585683 0.53471945 0.25368426 > > Due to this results I can't understand why...
Please note that optim() is (by default) *minimizing* the provided function, while maxLik() is *maximizing* the provided log-likelihood function. Cheers, Arne > A 26/09/2013, às 05:28, Arne Henningsen escreveu: > > Dear Filipe > > On 25 September 2013 14:23, Filipe Ribeiro <flipjribe...@hotmail.com> wrote: > > Hello everybody! > > > I'm having some trouble to compute maximum likelihood > > estimations using maxLik package and I hope that you > > could give me a hint. > > The main problem is that I'm not able to get a result not > > even close to the ones given by glm() directly, and the > > second one is: "Error in maxNRCompute(fn = logLikAttr, > > fnOrig = fn, gradOrig = grad, hessOrig = hess, : NA in > > gradient". > > > The codes: > > loglike.GGompiMaxLik <- function(theta,age,deaths,exposures) { > > alpha <- exp(theta[1]) > > beta <- exp(theta[2]) > > gamma <- exp(theta[3]) > > first <- alpha*exp(beta*age) > > second <- 1+(((alpha*gamma)/beta)*(exp(beta*age)-1)) > > mu <- first/second > > llk <- -sum((deaths * log(mu)) + (- mu*exposures)) > > return(llk) > > } > > > > fit1 <- maxLik(loglike.GGompiMaxLik, > > age=0:6, > > deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198), > > exposures=c(935023.67, 819392.00, 724568.17, 470947.00, > > 231951.64, 69502.65, 15798.72), > > start=c(-4.1402817, -0.6375773, -1.6945914)) > > > Do you know how I can solve this problem? > > > You did not write which model specification you want to estimate but I > am pretty sure that something in your log-likelihood function is > incorrect. The log-likelihood value at the starting values of the > parameters is so large that R even cannot calculate the likelihood > value: > > a <- loglike.GGompiMaxLik(c(-4.1402817, -0.6375773, -1.6945914), age=0:6, > > + deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198), > + exposures=c(935023.67, 819392.00, 724568.17, 470947.00, > + 231951.64, 69502.65, 15798.72)) > > a > > [1] 580365.2 > > exp(a) > > [1] Inf > > In the second iteration, the first parameter gets so small (large in > absolute terms, -5e+10) that the log-likelihood value become extremely > (numerically infinitely) large and the gradients cannot be computed > (by the finite-difference method): > > fit1 <- maxLik(loglike.GGompiMaxLik, > > + age=0:6, > + deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198), > + exposures=c(935023.67, 819392.00, 724568.17, 470947.00, > + 231951.64, 69502.65, 15798.72), > + start=c(-4.1402817, -0.6375773, -1.6945914)) > Iteration 2 > Parameter: > [1] -5.174233e+10 -3.839076e+02 5.988668e+00 > Gradient: > [,1] [,2] [,3] > [1,] NaN NaN NaN > Error in maxNRCompute(fn = logLikAttr, fnOrig = fn, gradOrig = grad, > hessOrig = hess, : > NA in gradient > > b <- loglike.GGompiMaxLik(c(-5.174233e+10, -3.839076e+02, 5.988668e+00), > age=0:6, > > + deaths=c(15545, 21278, 32444, 36201, 30360, 14201, 5198), > + exposures=c(935023.67, 819392.00, 724568.17, 470947.00, > + 231951.64, 69502.65, 15798.72)) > > b > > [1] Inf > > Please note that you can also find hints and ask questions about the > maxLik package in the forums at maxLik's R-Forge site: > > https://r-forge.r-project.org/projects/maxlik/ > > ...and please do not forget to cite the maxLik package in your publications: > > http://cran.r-project.org/web/packages/maxLik/citation.html > > Best wishes, > Arne -- Arne Henningsen http://www.arne-henningsen.name ______________________________________________ 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.