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
On 15-11-07 12:12 PM, Deniz OZONUR wrote:> 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-help at 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.
>