How to find the maximum likelihood estimators for this small data set?

[Vilter]

Let X_1,...,X_n be a random sample from a gamma(\alpha,\beta) population.

If \alpha and \beta are both unknown, there is no explicit formula for the MLEs of \alpha and \beta, but the maximum can be found numerically. Find the MLEs for \alpha and \beta for the below data.

First, I know how to find the MLE of \alpha and \beta. The rule is the sample average is the MLE of \alpha \beta. For each fixed value of \alpha, the value of \beta that maximizes L is \Sigma_i x_i/ (n \alpha). Substitute this into L. Then we just need to maximize the function of the one variable \alpha given by

\frac{1}{\Gamma(\alpha)^n (\Sigma_i x_i/ (n \alpha))^{n \alpha} } [\prod_{i=1}^n x_i]^{\alpha-1} e^{-n\alpha}

The data is:

22.0, 23.9, 20.9, 23.8, 25.0, 24.0, 21.7, 23.8, 22.8, 23.1, 23.1, 23.5, 23.0, 23.0

My attempt:

data <- c(22.0,23.9,20.9,23.8,25.0,24.0,21.7,23.8,22.8,23.1,23.1,23.5,23.0,23.0)

mean(data)

f = function(x){
  1/ ((gamma(x))^14 * (23.11429/x)^(14*x)) * prod(data)^(x-1) * exp(-14*x)
}

ans = optimize(f, interval = c(0,1000), maximum = TRUE)

# extract the coordinates of the maximum
x_max = ans$maximum
y_max = ans$objective

But my code gives me a wrong answer. I don't know where I got wrong. But mean(data) is correct.

[FJCC]

I would approach this in the following way. Be warned that I'm a novice at this sort of thing!
If the mean of the data is about 23 and is equal to \alpha\beta and the variance is about 1 and is equal to \alpha\beta^2, you can easily work out that \alpha is above 400. That results in gamma(400) returning Inf and other quantities in your f() function being very large. To get around that, I would work with the log-likelihood function. You can find that with a web search or just take the natural log of your calculation within f(). In calculating the log-likelihood, you will probably need to use the lgamma() function that returns log(gamma(x)).

[startz]

Adding to @FJCC's excellent point. If you have \alpha\beta=mean(data) and \alpha\beta^2=var(data), then the method of moments estimator is \beta=var(data)/mean(data) and \alpha=mean(data)/\beta.

These should give you good starting values. (And do be sure to follow @FJCC's advice about working with the log-likelihood.