2
$\begingroup$

Consider the following code in R:

data = c(603, 103, 225, 201, 1445, 1077, 1309, 6085, 469, 2309)
f = fitdistr(data, densfun="gamma")
rate = f$estimate["rate"]
loglik = f$loglik

This generates the following error:

Error in stats::optim(x = c(603, 103, 225, 201, 1445, 1077, 1309, 6085,  : 
  non-finite finite-difference value [1]
In addition: Warning messages:
1: In densfun(x, parm[1], parm[2], ...) : NaNs produced
2: In densfun(x, parm[1], parm[2], ...) : NaNs produced

Other posts here indicate that this is due to the extreme scaling of the data leading to unwanted 0's and infinities in the optimization routine:

https://stackoverflow.com/questions/15974773/errors-while-trying-to-fit-gamma-distribution-with-r-fitdistrmass

I note from these posts that the following workaround (scaling the data, then scaling back the rate) is possible:

f = fitdistr(data / 10, densfun="gamma")  ## scaled the data by 0.1
rate = f$estimate["rate"] / 10  ## scaled the rate by 0.1 - all good!

However it's vital to the future of mankind that I accurately calculate the log-likelihood of the fit too. Unfortunately the scaling of the data affects the loglik calculation and I cannot find a way to calculate it for the unscaled dataset:

loglik = f$loglik $$ THIS IS WRONG NOW!! :(

Can you suggest:

1) Parameters to fitdistr which might work around the problem?

2) An alternative distribution fitting library for R which might not suffer from the original problem?

3) A quick and easy alternative approach in a non-R environment to do the same job? (Python, Matlab etc)

|cite|improve this question
$\endgroup$
  • $\begingroup$ This works: library(help="MASS"); gdata = c(603, 103, 225, 201, 1445, 1077, 1309, 6085, 469, 2309); gfit=glm(gdata~1,family=Gamma(link="identity")); mean=fitted(gfit)[1]; gamma.shape(gfit); scale=mean/gamma.shape(gfit)$alpha; rate=1/scale ... but even starting at the solution, fitdistr won't be happy without scaling. $\endgroup$ – Glen_b -Reinstate Monica Aug 18 '16 at 10:15
  • $\begingroup$ The reason for this is the numeric approximation of the Hessian (to obtain an estimate of the variance-covariance matrix later on). This is much more likely to run into numeric problems if the two parameter estimates differ by several orders of magnitude as in the unscaled problem. If you switch off the hessian = TRUE option in the fitdistr code (it's not an exported argument) then fitdistr(data, "gamma") obtains the correct point estimates. Note that your glm()-based solution also avoids computing the full Hessian. $\endgroup$ – Achim Zeileis Aug 18 '16 at 18:27
4
$\begingroup$

I would definitely recommend scaling your data, probably even with 1000 and not just 10. I can't think of any reason why it would be necessary (let alone vital to the future of mankind) to compute the log-likelihood on the original scale rather than the transformed scale. Simply do all calculations on the transformed scale and you are fine.

Having said that it is easy to simply call sum(dgamma(..., log = TRUE)) with the transformed parameters to obtain the log-likelihood on the desired scale. For illustration let's look at the fit on the scale divided by 10 as you did and divided by 1000 as I would recommend:

## scaled data
d10 <- data/10
d1000 <- data/1000
## corresponding fitdistr
f10 <- fitdistr(d10, "gamma")
f1000 <- fitdistr(d1000, "gamma")

Then we can look at the shape parameters which essentially agree

s10 <- coef(f10)["shape"]
s1000 <- coef(f1000)["shape"]
c(s10, s1000)
##     shape     shape 
## 0.8971837 0.8973881

And then we compute the rate parameters on the scale divided by 10 and they also essentially agree:

r10_10 <- coef(f10)["rate"]
r10_1000 <- coef(f1000)["rate"]/100
c(r10_10, r10_1000)
##       rate       rate 
## 0.00648812 0.00649091

And then we can compute the log-likelihoods by hand or extract them from f10. 

l10_10 <- sum(dgamma(d10, shape = s10, rate = r10_10, log = TRUE))
l10_1000 <- sum(dgamma(d10, shape = s1000, rate = r10_1000, log = TRUE))
c(l10_10, l10_1000, logLik(f10))
## [1] -59.25158 -59.25158 -59.25158

They also essentially agree but the one from the fit on the scale divided by 1000 leads to a very slightly higher (= better) log-likelihood. This reflects that the numerical properties are better on this scale as both parameters (and their standard errors) are in the same order of magnitude. 

c(l10_10, logLik(f10)) < l10_1000
## [1] TRUE TRUE

To conclude we can also compute the log-likelihood on the original scale:

r1_10 <- coef(f10)["rate"]/10
r1_1000 <- coef(f1000)["rate"]/1000
l1_10 <- sum(dgamma(data, shape = s10, rate = r1_10, log = TRUE))
l1_1000 <- sum(dgamma(data, shape = s1000, rate = r1_1000, log = TRUE))
c(l1_10, l1_1000)
## [1] -82.27743 -82.27743

Of course, this differs from logLik(f10) and logLik(f1000) just by constants and all inference (e.g., likelihood ratio tests) or information criteria would lead to idential conclusions on any of the scales.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ Thanks for that, manually calculating the loglik using "sum(dgamma(..." is exactly what I needed, and I learnt a bit along the way too. Sorry did I say "future of mankind"? I meant "my assignment due next week". $\endgroup$ – Brendan Hill Aug 19 '16 at 2:35
  • $\begingroup$ It was clear that you were exaggerating but I still did not understand why it should be necessary to perform the task on the original scale. All questions such as inference or model selection can be answered based on the scaled data equally well. $\endgroup$ – Achim Zeileis Aug 19 '16 at 7:27
0
$\begingroup$

DIY (R code):

# - log likelihood
ll = function(par){
if(par[1]>0 & par[2]>0) return(-sum(dgamma(data,shape=par[1],rate=par[2],log=T)))
else return(Inf)
}

# initial point
init0 = c(1,0.1)

# Optimisation step
OPT = optim(init0,ll)

# -log likelihood
OPT$value

# MLE
OPT$par

# fitted density
fit.den = function(x) dgamma(x,shape=OPT$par[1],rate=OPT$par[2])

hist(data,probability=T)
curve(fit.den,0,8000,add=T,col="red")
|cite|improve this answer
$\endgroup$
  • 1
    $\begingroup$ This is very close to what fitdistr() does internally - which just tries to do more about choosing starting values and the method argument for the optim() call. The difference that causes the error is that your example does not compute the Hessian for the variance-covariance matrix. If you change your call to optim(init0, ll, hessian = TRUE) it also fails. The source is the scaling of the data that leads to a poorly conditioned Hessian. One could either use a better algorithm for the Hessian, e.g., numDeriv::hessian(ll, OPT$par), or simply scale the data appropriately (preferred). $\endgroup$ – Achim Zeileis Aug 18 '16 at 18:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.