Fitting a pareto distribution with (python) Scipy

Question

I have a data set that I know has a Pareto distribution. Can someone point me to how to fit this data set in Scipy? I got the below code to run but I have no idea what is being returned to me (a,b,c). Also, after obtaining a,b,c, how do I calculate the variance using them?

import scipy.stats as ss 
import scipy as sp

a,b,c=ss.pareto.fit(data)

Vergos

Answer 1

Be very careful fitting power laws!! Many reported power laws are actually badly fitted by a power law. See Clauset et al. for all the details (also on arxiv if you don't have access to the journal). They have a companion website to the article which now links to a Python implementation. Don't know if it uses Scipy because I used their R implementation when I last used it.

Nernst

The python implementation (code.google.com/p/agpy/wiki/PowerLaw) includes two versions; one depends on numpy, one does not. (I wrote it) – Denaturalize

and test fitting with this code to pareto and other distributions for chi-sqrt< 0.5 to reject H0, - the possibility that no association exists between the independent and dependent variables – Sigrid

Answer 2

Here's a quickly written version, taking some hints from the Reference page that Rupert gave. This is currently work in progress in scipy and statsmodels and requires MLE with some fixed or frozen parameters, which is only available in the trunk versions. No standard errors on the parameter estimates or other result statistics are available yet.

'''estimating pareto with 3 parameters (shape, loc, scale) with nested
minimization, MLE inside minimizing Kolmogorov-Smirnov statistic

running some examples looks good
Author: josef-pktd
'''

import numpy as np
from scipy import stats, optimize
#the following adds my frozen fit method to the distributions
#scipy trunk also has a fit method with some parameters fixed.
import scikits.statsmodels.sandbox.stats.distributions_patch

true = (0.5, 10, 1.)   # try different values
shape, loc, scale = true
rvs = stats.pareto.rvs(shape, loc=loc, scale=scale, size=1000)

rvsmin = rvs.min() #for starting value to fmin


def pareto_ks(loc, rvs):
    est = stats.pareto.fit_fr(rvs, 1., frozen=[np.nan, loc, np.nan])
    args = (est[0], loc, est[1])
    return stats.kstest(rvs,'pareto',args)[0]

locest = optimize.fmin(pareto_ks, rvsmin*0.7, (rvs,))
est = stats.pareto.fit_fr(rvs, 1., frozen=[np.nan, locest, np.nan])
args = (est[0], locest[0], est[1])
print 'estimate'
print args
print 'kstest'
print stats.kstest(rvs,'pareto',args)
print 'estimation error', args - np.array(true)

Elonore

Answer 3

Let's say you data is formated like this

import openturns as ot
data = [
    [2.7018013],
    [8.53280352],
    [1.15643882],
    [1.03359467],
    [1.53152735],
    [32.70434285],
    [12.60709624],
    [2.012235],
    [1.06747063],
    [1.41394096],
]
sample = ot.Sample([[v] for v in data])

You can easily fit a Pareto distribution using ParetoFactory of OpenTURNS library:

distribution = ot.ParetoFactory().build(sample)

You can of course print it:

print(distribution)
>>> Pareto(beta = 0.00317985, alpha=0.147365, gamma=1.0283)

or plot its PDF:

from openturns.viewer import View

pdf_graph = distribution.drawPDF()
pdf_graph.setTitle(str(distribution))
View(pdf_graph, add_legend=False)

More details on the ParetoFactory are provided in the documentation.

Franciscka

I tried running your code, it does not work for me. First of all, I think it should be distribution = ot.ParetoFactory.build(data). Secondly, I get a TypeError: Wrong number or type of arguments for overloaded function 'ParetoFactory_build'. It looks like Openturns library wants a specific format for 'sample', but I don't understand how it should be formatted. – Desman

Thank you @Desman for your feedback. It should be: distribution = ot.ParetoFactory().build(data). I corrected the code above. To be understood by OpenTURNS library, a Sample should be shaped (size, dimension). In my example it is (10, 1). – Franciscka

It seems like 'ot.ParetoFactory().build(sample)' struggles with zeros, returning "RuntimeError: InternalException : Ceres terminated with failure.". If I add a large enough number to the data with 'sample = ot.Sample(data + 1)' (one seems to work but a number less than one doesn't always work), the program runs with the number added to the gamma parameter. – Numerical

Answer 4

Before passing the data to build() function in OPENTURNS, make sure to convert it this way:

data = [[i] for i in data]

Because Sample() function may return an error.

FYI @Tropilio

Celestinecelestite

How does this answer the question? – Slump