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?

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)

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.