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@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.
