Source of the question: the Exercise 7.2(b) of Statistical Inference Book by Casella and Berger.
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.
A better way is to rather consider the log likelihood, which tends to simplify things. I urge you to read about likelihood estimation atleast to some degree to understand why this is customary. In this case is given by
$$ \ln(\mathcal{L(k,\theta}) = (k-1)\sum_{i=1}^N \ln(x_i)-\sum_{i=1}^N \frac{x_i}{\theta}-Nk\ln(\theta)-N\ln(\Gamma(k)) $$
as stated in wikipedia. Solving for $\frac{d}{dk}\ln(\mathcal{L(k,\theta})=0$ and $\frac{d}{d\theta}\ln(\mathcal{L(k,\theta})=0$ gives
$$ \hat{\theta}=\frac{1}{kN}\sum_{i=1}^N x_i $$
and
$$ \ln(\hat{k})-\psi(\hat{k})=\ln\left(\frac{1}{N}\sum_{i=1}^Nx_i \right) - \frac{1}{N}\sum_{i=1}^N \ln(x_i) $$
where the latter is solved numerically and $\psi(k)$ is the digamma function. If you plug in the numbers, you should get $\hat{k}\approx23.1$ and $\hat{\theta}\approx0.04$.