Statistical Distribution Fitter

This particular code is modified source from the pyeq2.  pyeq2 provides the backend that powers the ZunZun curve fitting site.  It's a fantastic site; a rather unique one which delivers well fitted curves for a given data set, fit statistics, and inspiring Bible quotes all together. Taking the code statistical distribution fitting code I was able to free this fragment, and made some improvements for debugging purposes.

Fitting distributions to data points is a common operation for me, because I commonly analyze data from different sources, apply wavelet filtering and so on.  Frequently I find myself in need of fitting a distribution to a given data set, I don't want to bog down ZunZun with requests.  Future improvements might include computing the K-S statistic along with the AIC criteria that is computed below.  Implementers may want to check the exec line and ensure that the scipy.statistics module is loaded under the alias spst.  Else the source line can be changed to match the requirements.

Statistical Distribution Fitter Code

# directly passing the distribution instance can yield "can't pickle instancemethod"
# exceptions, so the distribution name is passed instead
def SolveStatisticalDistribution(distributionName, data, inCriteriaForUseInListSorting):
    criteriaList = ['AIC', 'AICc_BA', 'nnlf']
    if inCriteriaForUseInListSorting not in criteriaList:
        raise Exception('Criteria to calculate for use in sorting was not in', str(criteriaList))

    try:
        exec('distribution = spst.' + distributionName)
    except:
        print 'Unable to find distribution, %s!' % distributionName
        return 0

    # only need to calculate these once
    eps = np.finfo(float).eps * 2.0
    data_min = data.min()
    data_max = data.max()
    data_mean = data.mean()
    data_range = data_max - data_min
    data_std_dev = np.std(data, dtype=np.float32) # note on precision http://docs.scipy.org/doc/np/reference/generated/np.std.html        
    
        
    # Try different starting parameters
    best_nnlf = 1.0E300
    best_parameters = None

    if distribution.name in ['beta']:
        try:
            rangeData = (data - data_min) / data_range
            data_mean = rangeData.mean()
            data_var = rangeData.var()
            
            par_a = data_mean * ((data_mean * (1.0 - data_mean) / data_var) - 1.0)
            par_b = (1.0 - data_mean) * ((data_mean * (1.0 - data_mean) / data_var) - 1.0)
            par0 = (par_a, par_b, data_min - 0.001*data_min, data_range * 1.001)
            par_est = tuple(distribution.fit(data, *(par0[:-2]), loc = data_min - 0.01*data_min, scale = data_range * 1.01))
            nnlf = distribution.nnlf(par_est, data)
            if np.isfinite(nnlf) and nnlf < best_nnlf and nnlf > 0.0:
                best_parameters = par_est
                best_nnlf = nnlf
        except:
            print 'Failure fitting a beta!'

    if distribution.name in ['truncnorm','betaprime','reciprocal']:
        try:
            par0 = (data_mean-2.0*data_std_dev, data_mean+2.0*data_std_dev)
            par_est = tuple(distribution.fit(data, loc=data_mean, scale=data_std_dev, *par0))
            nnlf = distribution.nnlf(par_est, data)
            if np.isfinite(nnlf) and nnlf < best_nnlf and nnlf > 0.0:
                best_parameters = par_est
                best_nnlf = nnlf
        except:
            print 'Failure fitting a trunk norm!'
        
    try:
        par_est = tuple(distribution.fit(data))
        nnlf = distribution.nnlf(par_est, data)
        if np.isfinite(nnlf) and nnlf < best_nnlf and nnlf > 0.0:
            best_parameters = par_est
            best_nnlf = nnlf
    except:
        print 'Failure fitting no values!'

    try:
        par_est = tuple(distribution.fit(data, loc=0.0, scale=1.0))
        nnlf = distribution.nnlf(par_est, data)
        if np.isfinite(nnlf) and nnlf < best_nnlf and nnlf > 0.0:
            best_parameters = par_est
            best_nnlf = nnlf
    except:
        print 'Failure to fit location and scale at default!'
        
    try:
        par_est = tuple(distribution.fit(data, loc=data_mean, scale=data_std_dev))
        nnlf = distribution.nnlf(par_est, data)
        if np.isfinite(nnlf) and nnlf < best_nnlf and nnlf > 0.0:
            best_parameters = par_est
            best_nnlf = nnlf
    except:
        print 'Failure fitting with mean as location and std as scale!'

    try:
        par_est = tuple(distribution.fit(data, loc=data_max+eps, scale=data_std_dev))
        nnlf = distribution.nnlf(par_est, data)
        if np.isfinite(nnlf) and nnlf < best_nnlf and nnlf > 0.0:
            best_parameters = par_est
            best_nnlf = nnlf
    except:
        print 'Failure fitting a +max eps over std!'

    try:
        par_est = tuple(distribution.fit(data, loc=data_min-eps, scale=data_std_dev))
        nnlf = distribution.nnlf(par_est, data)
        if np.isfinite(nnlf) and nnlf < best_nnlf and nnlf > 0.0:
            best_parameters = par_est
            best_nnlf = nnlf
    except:
        print 'Failure fitting a min-eps over std!'

    try:
        par_est = tuple(distribution.fit(data, loc=data_max+eps, scale=data_range))
        nnlf = distribution.nnlf(par_est, data)
        if np.isfinite(nnlf) and nnlf < best_nnlf and nnlf > 0.0:
            best_parameters = par_est
            best_nnlf = nnlf
    except:
        print 'Failure fitting a +eps over range!'

    try:
        par_est = tuple(distribution.fit(data, loc=data_min-eps, scale=data_range))
        nnlf = distribution.nnlf(par_est, data)
        if np.isfinite(nnlf) and nnlf < best_nnlf and nnlf > 0.0:
            best_parameters = par_est
            best_nnlf = nnlf
    except:
        print 'Failure fitting a -min-eps over range!'

    if (best_nnlf < 1.0E300) and (best_parameters is not None):
        try:
            k = len(best_parameters)
            AIC = 2.0*k + 2.0 * best_nnlf
            n = len(data)
            AICc_BA = AIC + ( (2.0 * k * (k+1.0)) / (n - k - 1.0))

            temp = {}
            temp['distributionName'] = distributionName
            temp['fittedParameters'] = best_parameters
            temp['nnlf'] = best_nnlf
            temp['AIC'] = AIC
            temp['AICc_BA'] = AICc_BA
            
            if inCriteriaForUseInListSorting == 'nnlf':
                return [best_nnlf, temp]
            elif inCriteriaForUseInListSorting == 'AIC':
                return (AIC, temp)
            else:
                return (AICc_BA, temp)
        except:
            print 'Exception in final calculations!'
    else:
        print 'Failed to find best parameters!'
        return 0