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.