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:
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)
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:15hessian = TRUE
option in thefitdistr
code (it's not an exported argument) thenfitdistr(data, "gamma")
obtains the correct point estimates. Note that yourglm()
-based solution also avoids computing the full Hessian. $\endgroup$ – Achim Zeileis Aug 18 '16 at 18:27