This question might also be suited to the programming site, but I thought since there is enough on the statistics side, I would use this forum. I am trying to estimate the alpha parameter in a Gamma distribution using maximum likelihood method, and using the optimization functions available in R.
To begin with, I generated a random sample from Gamma(Alpha, Beta) in R.
shape <- 2
scale <- 1.5
set.seed(123456)
myData <- round(rgamma(n=50, shape=shape, scale=scale),2)
Using the maximum likelihood estimation method, and setting up the likelihood function to be in terms of alpha only, I created a function in R and I am trying to optimize it. So I wrote the likelihood function, took the log, took the partial derivative with respect to Beta, and found the MLE of Beta. I then substituted the MLE of Beta back into the likelihood function to arrive at the likelihood in terms of alpha only. My function is as follows:
objFunction <- function(myData, alpha) {
sumX <- sum(myData)
prodX <- prod(myData)
n <- length(myData)
estimate <- (1/((gamma(alpha^n))*((sumX/(n*alpha))^(n*alpha))))*((prodX)^(alpha-1))*(exp(1)^(-n*alpha))
return(-1*estimate)
}
Now to optimize, I attempted three different functions from R:
optim(par=0, fn=objFunction, method = "Brent", lower = 0, upper = 10, alpha=2)
nlm(objFunction, momAlpha, myData=myData)
optimize(f=objFunction, c(0,10), alpha=2, maximum=TRUE)
The variable momAlpha, is basically the method of moments estimator for the Alpha, as that would be a good start. Just for completeness:
momAlpha <- (mean(myData)^2)/var(myData)
momBeta <- var(myData)/mean(myData)
These are available in many online references.
Now when I ran the optimization functions above, my results were not clear to me and I need some help understanding:
optim(par=0, fn=objFunction, method = "Brent", lower = 0, upper = 10, alpha=2)
$par
[1] 0.000000005349424
$value
[1] -101196146
$counts
function gradient
NA NA
$convergence
[1] 0
$message
NULL
Why is this estimate way out of range?
nlm(objFunction, momAlpha, myData=myData)
$minimum
[1] 0
$estimate
[1] 1.919078
$gradient
[1] 0
$code
[1] 1
$iterations
[1] 0
Warning messages:
1: In f(x, ...) : value out of range in 'gammafn'
2: In f(x, ...) : value out of range in 'gammafn'
3: In f(x, ...) : value out of range in 'gammafn'
The estimate here is nothing but the starting point I provided, why?
optimize(f=objFunction, c(0,2), alpha=2, maximum=TRUE)
$maximum
[1] 1.999934
$objective
[1] -0.2706795
Is this even right?
I am still developing my intuition for the subject, but it seems that I am either doing something wrong, i.e. my objective function is incorrect or the parameter settings of the functions is incorrect or I simply don't understand the way the functions work. I appreciate any help in guiding through this!