[R] fitdistr, mle's and gamma distribution

Spencer Graves spencer.graves at pdf.com
Tue Sep 30 22:07:11 CEST 2003


      In my experience, the most likely cause of this problem is that 
optim may try to test nonpositive values for shape or scale.  I avoid 
this situation by programming the log(likelihood) in terms of log(shape) 
and log(scale) as follows: 

 > gammaLoglik <-
+ function(x, logShape, logScale, negative=TRUE){
+ lglk <- sum(dgamma(x, shape=exp(logShape), scale=exp(logScale),
+ log=TRUE))
+ if(negative) return(-lglk) else return(lglk)
+ }
 > tst <- rgamma(10, 1)
 > gammaLoglik(tst, 0, 0)
[1] 12.29849

Then I then call optim directly to minimize the negative of the 
log(likelihood). 

      If I've guessed correctly, this should fix the problem.  If not, 
let us know. 

      hope this helps.  spencer graves

Lourens Olivier Walters wrote:

>Dear R Users, 
>
>I am trying to obtain a best-fit analytic distribution for a dataset
>with 11535459 entries. The data range in value from 1 to 300000000. I
>use: fitdistr(data, "gamma") to obtain mle's for the parameters.
> 
>I get the following error:
>
>Error in optim(start, mylogfn, x = x, hessian = TRUE, ...) :
>        non-finite finite-difference value [1]
>
>And the following warnings:
>
>NaNs produced in: dgamma(x, shape, scale, log)
>
>I have the same problem with the exponential distribution, but the
>lognormal and weibull distributions don't have this problem. 
>
>I suspect the problem has to do with:
>
>"It means that in computing derivatives by finite differencing, one 
>of the values is NA, +Inf or -Inf. Put some print values in 
>your objective function and find out why it is giving a non-finite
>value." - quote by Prof Ripley from r-help
>
>as from the warnings it is clear that some of the values obtained during
>maximum likelihood computations are NaN. I cannot however print values
>of the objective function, as in my case it is dgamma which is called by
>fitdistr, over which I don't have control.  
>
>Can anyone please point me in the right direction?
>
>Lourens
>
>  
>




More information about the R-help mailing list