I'm trying to calculate the two parameters of the Gamma distribution by solving the system of two equations obtained by differentiating the Log-likelihood function relatively to both parameters as suggested here Maximum likelihood estimators for Gamma distribution.
Here is the PDF (Probability Density Function) of the Gamma distribution I'm using (similarly to the one given in the link mentioned above):
$$f(x) = \frac {1}{\Gamma(r)}\lambda^{r}x^{r-1}e^{-\lambda x} \quad \text{with} \quad x \geq 0 \quad \text{and} \quad r, \lambda > 0 $$
where:
- $r$ is the shape (= $\alpha$).
- $\lambda$ is the rate (= $\frac{1}{\beta}$ where $\beta$ is the scale).
And, the Log-likelihood to maximize is given by:
$$logL=(r-1)\sum_ilogx_i-\lambda T +(nr)log\lambda -nlog(\Gamma(r))$$
where:
- $T= x_1 + x_2 + ... + x_n $
Then I got the following system of two equations from the derivatives of $logL$, for which I'm trying to calculate the roots:
$$ \frac{\partial logL}{\partial \lambda} = -T + \frac{nr}{\lambda} = 0 $$ $$ \frac{\partial logL}{\partial r} = \sum_ilogx_i + nlog\lambda - n\frac{\partial \log(\Gamma(r))}{\partial r} = 0 $$
Now knowing that $\frac{\partial \log(\Gamma(r))}{\partial r}$ is equal to the Digamma function which can be approximated (according to Sympy) with the Euler-Mascheroni constant, how can I solve the previous system of two equations in order to calculate $(r, \lambda)$?
If you are approximating that term by a constant, solve the second equation for $\lambda$ and then use it in the first to get $r$.
The approximation removes $r$ from the second equation.
However, I think a constant isn't a very good approximation, is it?