Re: [R] MLE Estimation of Gamma Distribution Parameters for data with 'zeros'

From: Peter Dalgaard <p.dalgaard_at_biostat.ku.dk>
Date: Wed, 11 Jun 2008 22:40:58 +0200

Fox, Aaron wrote:
> Greetings, all
>
> I am having difficulty getting the fitdistr() function to return without
> an error on my data. Specifically, what I'm trying to do is get a
> parameter estimation for fracture intensity data in a well / borehole.
> Lower bound is 0 (no fractures in the selected data interval), and upper
> bound is ~ 10 - 50, depending on what scale you are conducting the
> analysis on.
>

> I read in the data from a text file, convert it to numerics, and then
> calculate initial estimates of the shape and scale parameters for the
> gamma distribution from moments. I then feed this back into the
> fitdistr() function.
>
> R code (to this point):
> ########################################
> data.raw=c(readLines("FSM_C_9m_ENE.inp"))
> data.num <- as.numeric(data.raw)
> data.num
> library(MASS)
> shape.mom = ((mean(data.num))/ (sd(data.num))^2
> shape.mom
> med.data = mean(data.num)
> sd.data = sd(data.num)
> med.data
> sd.data
> shape.mom = (med.data/sd.data)^2
> shape.mom
> scale.mom = (sd.data^2)/med.data
> scale.mom
> fitdistr(data.num,"gamma",list(shape=shape.mom,
> scale=scale.mom),lower=0)
>
> fitdistr() returns the following error:
>
> " Error in optim(x = c(0.402707037, 0.403483333, 0.404383704,
> 2.432626667, :
> L-BFGS-B needs finite values of 'fn'"
>
> Next thing I tried was to manually specify the negative log-likelihood
> function and pass it straight to mle() (the method specified in Ricci's
> tutorial on fitting distributions with R). Basically, I got the same
> result as using fitdistr().
>
> Finally I tried using some R code I found from someone with a similar
> problem back in 2003 from the archives of this mailing list:
>
> R code
> ########################################
> gamma.param1 <- shape.mom
> gamma.param2 <- scale.mom
> log.gamma.param1 <- log(gamma.param1)
> log.gamma.param2 <- log(gamma.param2)
>
>
> gammaLoglik <-
> function(params,
> negative=TRUE){
> lglk <- sum(dgamma(data, shape=exp(params[1]), scale=exp(params[2]),
> log=TRUE))
> if(negative)
> return(-lglk)
> else
> return(lglk)
> }
>
> optim.list <- optim(c(log.gamma.param1, log.gamma.param2), gammaLoglik)
> gamma.param1 <- exp(optim.list$par[1])
> gamma.param2 <- exp(optim.list$par[2])
> #########################################
>
> If I test this function using my sample data and the estimates of shape
> and scale derived from the method of moments, gammaLogLike returns as
> INF. I suspect the problem is that the zeros in the data are causing the
> optim solver problems when it attempts to minimize the negative
> log-likelihood function.
>
> Can anyone suggest some advice on a work-around? I have seen
> suggestions online that a 'censoring' algorithm can allow one to use MLE
> methods to estimate the gamma distribution for data with zero values
> (Wilkes, 1990, Journal of Climate). I have not, however, found R code to
> implement this, and, frankly, am not smart enough to do it myself... :-)
>
> Any suggestions? Has anyone else run up against this and written code to
> solve the problem?
>
It's fairly easy. You decide that the zeros really represent values less than "delta" (e.g. 0.5 if your data are integers), then replace dgamma(0,....) with pgamma(delta,...) in the likelihood. (And, BTW, the problem is not that the optimizers get in trouble, but rather that the log-likelihood *is* +/- Inf if there are zeros in data unless the shape parameter is exactly 1 -- the x^(a-1) factor in the gamma density causes this).

-- 
   O__  ---- Peter Dalgaard             ุster Farimagsgade 5, Entr.B
  c/ /'_ --- Dept. of Biostatistics     PO Box 2099, 1014 Cph. K
 (*) \(*) -- University of Copenhagen   Denmark      Ph:  (+45) 35327918
~~~~~~~~~~ - (p.dalgaard_at_biostat.ku.dk)              FAX: (+45) 35327907

______________________________________________
R-help_at_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.
Received on Wed 11 Jun 2008 - 21:04:25 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Wed 11 Jun 2008 - 21:31:07 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive