使用（python）Scipy拟合帕累托分布

Question

我有一个数据集，我知道它符合帕累托分布。有人能告诉我怎么在Scipy中拟合这个数据集吗？我运行了下面的代码，但我不知道返回的(a,b,c)是什么意思。另外，得到a、b、c之后，我该怎么用它们来计算方差呢？

import scipy.stats as ss 
import scipy as sp

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

Answer 1

假设你的数据格式是这样的：

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])

你可以很简单地使用OpenTURNS库中的ParetoFactory来拟合一个帕累托分布：

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

当然，你也可以打印出来：

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

或者画出它的概率密度函数（PDF）：

from openturns.viewer import View

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

关于ParetoFactory的更多细节，可以查看文档中的说明：ParetoFactory。

Answer 2

这是一个快速写的版本，参考了Rupert提供的参考页面的一些提示。 目前这个功能还在scipy和statsmodels中开发中，需要使用最大似然估计（MLE）并且有一些参数是固定或冻结的，这个功能目前只在开发版本中可用。 目前还没有关于参数估计的标准误差或其他结果统计的数据。

'''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)

Answer 3

在拟合幂律时要非常小心！！ 许多报告的幂律其实并没有被很好地拟合。想了解更多细节，可以查看Clauset等人的研究（如果你无法访问期刊，也可以在arxiv上找到）。他们还有一个配套网站，里面有与文章相关的内容，现在还链接到了一个Python的实现。我不确定它是否使用了Scipy，因为我上次用的是他们的R实现。