Using python scipy to fit gamma distribution to data

Question

I want to fit a gamma distribution to my data, which I do using this

import scipy.stats as ss
import scipy as sp
import numpy as np
import os
import matplotlib.pyplot as plt

alpha = []
beta = []
loc = []

data = np.loadtxt(data)
fit_alpha, fit_loc, fit_beta = ss.gamma.fit(data, floc=0, fscale=1)

I want to keep one of the parameters to the gamma distribution as a variable (say the shape), and fix one of the parameters (say scale=1). However, if I keep the loc variable as zero, I am not able to fix the scale at one. Is there some workaround for this? Can I not parametrize the gamma distribution using only the shape and scale?

Autocatalysis

You have run into a bug that has been fixed (github.com/scipy/scipy/pull/2519) in the development version of scipy. The fix will be in version 0.13 when it is released. – Maryjomaryl

Answer 1

In a comment I said you have run into a bug in the gamma distribution--it does not let you fix both the location and the scale. The bug was fixed in scipy 0.13, but if you can't upgrade, you can work around the bug by using the fit method of the class rv_continuous, which is the parent class of gamma:

In [22]: from scipy.stats import rv_continuous, gamma

In [23]: x = gamma.rvs(2.5, loc=0, scale=4, size=1000)  # A test sample.

In [24]: rv_continuous.fit(gamma, x, floc=0, fscale=4)
Out[24]: (2.5335837650122608, 0, 4)

Maryjomaryl

+1 for floc, documentation at docs.scipy.org/doc/scipy/reference/generated/… is devoid of any reference to this param. – Frodeen

Answer 2

Looking at the implementation of gamma.fit:

def fit(self, data, *args, **kwds):
    floc = kwds.get('floc', None)
    if floc == 0:
        xbar = ravel(data).mean()
        logx_bar = ravel(log(data)).mean()
        s = log(xbar) - logx_bar
        def func(a):
            return log(a) - special.digamma(a) - s
        aest = (3-s + math.sqrt((s-3)**2 + 24*s)) / (12*s)
        xa = aest*(1-0.4)
        xb = aest*(1+0.4)
        a = optimize.brentq(func, xa, xb, disp=0)
        scale = xbar / a
        return a, floc, scale
    else:
        return super(gamma_gen, self).fit(data, *args, **kwds)

If you put floc=None, it will call the fit function of the parent class (which is rv_continuous) and you can fix the scale.

Tannatannage