Fitting Gamma Parameters via MLE

Basic Concepts

We show how to estimate the parameters of the gamma distribution using the maximum likelihood approach. The pdf of the gamma distribution is

Gamma distribution pdf and soLikelihood function gamma distribution

It turns out that the maximum of L(α, β)  occurs when β = xÌ„ / α. We can now use Excel’s Solver to find the value of α that maximizes LL. Alternatively, we can use the following iteration method to find α.

Initializing the iteration

Iteration for gamma distribution

where ψ(z) (also denoted ψ0(z)) is the digamma function and ψ1(z) is the trigamma function. These functions can be estimated as follows:

Digamma function estimation

Trigamma function estimation

The estimates are quite accurate for values of z ≥ 4. For smaller values of z, we can improve the estimates by using the following properties:

Digamma trigamma function properties

Worksheet Function

Real Statistics Function: The digamma and trigamma functions can be computed via the following Real Statistic worksheet function:

POLYGAMMA(z, k) = digamma function at z if k = 0 (default) and trigamma function at z if k = 1.

Example

Example 1: Find the parameters of the gamma distribution that best fit the data in range A4:A18  of Figure 1.

The preliminary calculations are shown in range D4:D7 of Figure 1. Note that the formula in cell D7 is an array function (and so you must press Ctrl-Shft-Enter and not just Enter, unless you are using Excel 365). The iteration is shown in the range D14:D17. As you can see, the iteration converges quite rapidly.

Fitting gamma parameters MLE

Figure 1 – Fitting a Gamma Distribution

The alpha and beta parameters are 3.425 (cell D9) and 0.975 (cell D10).

Examples Workbook

Click here to download the Excel workbook with the examples described on this webpage.

References

Wikipedia (2017) Digamma function
https://en.wikipedia.org/wiki/Digamma_function

Wikipedia (2017) Trigamma function
https://en.wikipedia.org/wiki/Trigamma_function

12 thoughts on “Fitting Gamma Parameters via MLE”

  1. here is the pseudo code. thank u mr charles

    function psi(z,k)

    if z>=4 then
    if k=0 then
    psi = //digamma series
    exit function
    elseif k=1 then
    psi = //trigamma
    ` endif
    else
    psi = psi(z+1,k) – 1/z
    end if
    end function

    Reply
  2. Hey Charles!
    Thanks for your valuable contribution. Can you provide me with the codes to create the digamma and trigamma functions? I didnt understudy how to make them.
    Thanks
    Antonio

    Reply
  3. Hi charles,

    What formula do I have to put in my VBA? I don’t get it how ‘polygamma’ can be determined.

    Hope you can help me.

    Gr. Jonathan

    Reply
    • Hi Jonathan,
      The formulas for polygamma when k = 1 (digamma) and k = 2 (trigamma) are shown on this webpage. These are infinite sums, but you only need to use asmall number of the terms.
      Charles

      Reply
  4. This is great!
    If I understand correctly, you provided a method to estimate the shape and scale parameters for the gamma distribution. I was wondering, how would we estimate the rate (1/scale) parameter instead if that was our data?

    Thanks for your help.

    Reply
  5. Very helpful instruction that was easily implemented in Excel. THANK-YOU!
    Two comments:
    1. The nested formula for “mean ln x” imbeds a range in the ln function. This should be noted as an array formula (enter with CTRL SHIFT ENTER). If a simple ENTER is used, then the calculation of “mean ln x” is incorrect.
    2. POLYGAMMA function can be created as a user-defined function (UDF) using VBA. The given formulas (series expansions for z>4 and recursive relations for z<4) helped. But I had to look up the values for trigamma and digamma for 0<z<=1 and include these in my UDF for interpolating values in that range. This was critical to do correctly because the recursion relations start with decimal values of z if z is not an integer.

    Reply
    • Cal,
      I am pleased that the website has been helpful to your implementation.
      1. Thanks for identifying this omission. I have added this information to the webpage.
      2. I didn’t think that it was necessary to lookup values of digamma and trigamma for 0 < z <= 1. E.g. to calculate the digamma value at z = .2, you should be able to calculate the value psi(4.2) and then use the formula psi(z) = psi(z+1) - 1/z several times to get the value of psi(.2). psi(4.2) can be calculated as ln(4.2) - 1/(2*4.2) - 1/(12*4.2^2) + etc. psi(.2) = psi(1.2) - 1/.2 and psi(1.2) = psi(2.2) - 1/1.2 and psi(2.2) = psi(3.2) - 1/2.2 and psi(3.2) = psi(4.2) - 1/3.2. Thus once psi(4.2) is calculated, you should be able to calculate psi(.2). Thus z doesn't have to be an integer. Charles

      Reply

Leave a Comment