[R] About error: L-BFGS-B needs finite values of 'fn'

[Deniz OZONUR]

Hi,

I am trying to obtain power of Likelihood ratio test for comparing gamma distribution against generalized gamma distribution. And so I need maximum likelihood estimates of Generalized gamma distribution with three parameters. I wrote code as follows.

require(bbmle)
library("bbmle")

require(flexsurv)
library("flexsurv")

sig=0.05
den=1000
n=30
apar=2 ###alpha
bpar=3 ##beta
cpar=2 ##c parameter


LRatio=function(den,n,par=c(cpar,apar,bpar)){

LR2=rep(0,den)

count=rep(0,den)

cpar=par[1]
apar=par[2]
bpar=par[3]

for(i in 1:den){

y=rgengamma.orig(n,shape=cpar,scale=bpar,k=apar)

gamma4 = function(shape, scale) {
-sum(dgamma(y, shape = shape, scale = scale,log = TRUE))
}

gm = mean(y)
cv = var(y)/mean(y)

m5 = mle2(gamma4, start = list(shape = gm/cv, scale = cv),method = "L-BFGS-B", lower =c(.00001,.00001),upper = c(Inf,Inf))


gengamma3 = function(shape, scale,k) {
-sum(dgengamma.orig(y, shape = shape, scale = scale,k=k,log =TRUE))
}

ci=mean(y) #c initial value
a1=ci*mean(y)^(ci-1)
a2=ci*(ci-1)*(mean(y)^(ci-1))/2
mu1=mean(y)^ci+a2*mean(y^2)
mu2=(a1^2)*mean(y^2)+2*a1*a2*mean(y^3)+(a2^2)*(mean(y^4)-mean(y^2)^2)
alp =(mu1^2)/mu2 #alpha initial value
bet=mean(y)*gamma(alp)/gamma(alp+(1/ci)) #beta initial value

m6 = mle2(gengamma3, start = list(shape = ci, scale = bet, k=alp),method = "L-BFGS-B", lower = c(.00001,.00001,.00001),upper = c(Inf, Inf, nf))


LR2[i]=2*(logLik(m6)-logLik(m5))
count[i]=LR2[i]>=qchisq(1-sig, df=1)

}

pow=sum(count)/den
print(i)
print(pow)
}

But I got error : optim(par = c(3.88907163215354, 3.62005456122935, 1.66499331462506 : L-BFGS-B needs finite values of 'fn'


What is wrong? Can you hep me, thanks..
Deniz...

______________________________________________
R-h...@r-project.org mailing list -- To UNSUBSCRIBE and more, see
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.

[ProfJCNash]

Numerical gradient approximations are being used in your call, so my
guess is that the "epsilon" has made (parameter + epsilon) an
inadmissible argument for your likelihood. If you can supply analytical
gradients, the issue has a good chance of going away. Otherwise, you'll
need to use bounds or transformations to avoid the parameter region
giving undefined results.

JN

[Ravi Varadhan]

With a small sample size, n=30, you will have realizations of data where you will run into difficulties with the MLE of generalized Gamma distribution. This is mainly due to the `k' parameter. Increase the sample size (e.g., n=50 or 100) and this problem is less likely to happen (but can still happen).

I would strongly suggest that when you are doing simulations, you should encapsulate the parameter estimation inside a `try' or `tryCatch' statement so that when there is an error, the simulation keeps going rather than crashing out.

See the attached code.

Best,
Ravi

Ravi Varadhan, Ph.D. (Biostatistics), Ph.D. (Environmental Engg)
Associate Professor, Department of Oncology
Division of Biostatistics & Bionformatics
Sidney Kimmel Comprehensive Cancer Center
Johns Hopkins University
550 N. Broadway, Suite 1111-E
Baltimore, MD 21205
410-502-2619

[Ravi Varadhan]

It seems like there is substantial finite-sample bias in the MLEs. Either that or there is some error in your procedure. See attached code.

Ravi