-1

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)$?

Community
  • 1
Hakim
  • 59
  • Shouldn't that term be $-n\Gamma'(r)/\Gamma(r)$? – MPW Dec 26 '16 at 17:53
  • @MPW You're right, I've just corrected the formula. – Hakim Dec 26 '16 at 19:05
  • Looking at the graph of the Digamma function, I suspect its slope may be approximately constant for positive $x>>1$, but not the function itself – MPW Dec 26 '16 at 21:26
  • @MPW Indeed the derivative of the Digamma function is almost constant for large positive numbers, but does this help to solve the previous equations? – Hakim Dec 27 '16 at 15:49
  • 1
    If so, then you can approximate the Digamma function by $a+bx$ for some choice of $a$ and $b$. I think you are telling me the slope $b$ is the Euler-Mascheroni constant $\gamma$, not sure what to use for $a$. But just pick a couple of values and make a linear interpolation. – MPW Dec 27 '16 at 15:54
  • @MPW Actually, what I said earlier is that for a number $x$: $\psi(x) = -EulerMascheroniConstant + ax$ (according to sympy). This is due to the recurrence relation: $\psi(x+1) = \psi(x) + \frac{1}{x}$. – Hakim Dec 27 '16 at 18:15
  • @MPW Good news, I just noticed that the $f(x) = log(x)$ fits almost perfectly well with samples of the digamma function for $0 < x < 100$. – Hakim Dec 27 '16 at 18:29
  • @MPW The bad news is that if I replace the value of $r$ from the first equation and use $log(x)$ as an approximation for the Digamma function, I can't obtain a solution for $\lambda$ in the second equation. – Hakim Dec 27 '16 at 18:40

1 Answers1

0

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?

MPW
  • 43,638