[R] maximizing the gamma likelihood

markleeds at verizon.net markleeds at verizon.net
Fri May 23 23:28:25 CEST 2008


for learning purposes and also to help someone, i used roger peng's 
document to get the mle's of the gamma where the gamma is defined as

f(y_i) = (1/gammafunction(shape)) * (scale^shape) * (y_i^(shape-1)) * 
exp(-scale*y_i)

( i'm defining the scale as lambda rather than 1/lambda. various books 
define it differently  ).

i found the likelihood to be n*shape*log(scale)  + 
(shape-1)*sum(log(y_i) - scale*sum(y_i)
then i wrote below which is just roger peng's likelihood example but 
using the gamma instead of the normal.
I get estimates back but i separately found that the analytical mle of 
the scale is equal to 1/ analytical mle(shape).
and my estimates aren't consistent with that fact ? this leads me to 
assume that my estimates are not correct.

can anyone tell me what i'm doing wrong. maybe my starting values are 
too far off ? thanks.

make.negloglik <- function(data, fixed=c(FALSE,FALSE)) {
   op <- fixed
   function(p) {
     op[!fixed] <- p
     shape <- exp(op[1])
     scale <- exp(op[2])
     a <- length(data)*shape*log(scale)
     b <- (shape-1)*sum(log(data))
     c <- -1.0*scale*sum(data)
     -(a + b + c)
    }
}

vsamples<- c(14.7, 18.8, 14, 15.9, 9.7, 12.8)
nLL <- make.negloglik(vsamples)
temp <- optim(c(scale=1,shape=1), nLL, method="BFGS")[["par"]]
estimates <- log(temp)
print(estimates)

check <- estimates[1]/mean(vsamples)
print(check)



More information about the R-help mailing list