python scipy stats 帕累托拟合：它是如何工作的

Question

...帮助和在线文档说函数 scipy.stats.pareto.fit 将要拟合的数据集作为变量，以及可选的 b（指数）、loc、scale。 结果以三元组形式出现（指数、loc、标度）

从相同的分布生成数据应该导致拟合找到用于生成数据的参数，例如（使用python 3 colsole）

$  python
Python 3.3.0 (default, Dec 12 2012, 07:43:02) 
[GCC 4.7.2] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>>

（在下面的代码行中省略了 python 控制台提示符“>>>”）

dataset=scipy.stats.pareto.rvs(1.5,size=10000)  #generating data
scipy.stats.pareto.fit(dataset)

然而这会导致

(1.0, nan, 0.0)

（指数 1，应为 1.5）和

dataset=scipy.stats.pareto.rvs(1.1,size=10000)  #generating data
scipy.stats.pareto.fit(dataset)

结果

(1.0, nan, 0.0)

（指数 1，应为 1.1）和

dataset=scipy.stats.pareto.rvs(4,loc=2.0,scale=0.4,size=10000)    #generating data
scipy.stats.pareto.fit(dataset)

（指数应为 4，loc 应为 2，小数位数应为 0.4）in

(1.0, nan, 0.0)

等等。 调用 fit 函数时给出另一个指数

scipy.stats.pareto.fit(dataset,1.4)

始终返回该指数

(1.3999999999999999, nan, 0.0)

明显的问题是：我是否完全误解了这个拟合函数的目的，它的使用方式是否有所不同，或者只是被破坏了？

备注：之前有人提到像 Aaron Clauset 网页上给出的专用函数 (http://tuvalu.santafe.edu/~aaronc/powerlaws/) 比 scipy.stats 方法更可靠，应该使用相反：这可能是真的，但它们也非常非常非常非常耗时，并且在普通 PC 上处理 10000 个点的数据集需要花费很多个小时（可能是几天、几周、几年）。

编辑：哦：拟合函数的参数不是分布的指数而是指数负1（但这并不能改变上面的问题）

Answer 1

看来您必须对

loc

和

scale

提供猜测：

In [78]: import scipy.stats as stats

In [79]: b, loc, scale = 1.5, 0, 1

In [80]: data = stats.pareto.rvs(b, size=10000)

In [81]: stats.pareto.fit(data, 1, loc=0, scale=1)
Out[81]: (1.5237427002368424, -2.8457847787917788e-05, 1.0000329980475393)

并且猜测必须非常准确才能成功拟合：

In [82]: stats.pareto.fit(data, 1, loc=0, scale=1.01)
Out[82]: (1.5254113096223709, -0.0015898489208676779, 1.0015943893384001)

In [83]: stats.pareto.fit(data, 1, loc=0, scale=1.05)
Out[83]: (1.5234726749064218, 0.00025804526532994751, 0.99974649559141171)

In [84]: stats.pareto.fit(data, 1, loc=0.05, scale=1.05)
Out[84]: (1.0, 0.050000000000000003, 1.05)

希望问题的背景能够告诉您对

loc

和

scale

的适当猜测应该是什么。最有可能的是，

loc=0

和

scale=1

。

Answer 2

拟合方法是一种非常通用且简单的方法，它对分布的非负似然函数（self.nnlf）进行optimize.fmin。 在像帕累托这样的分布中，其参数可以创建未定义的区域，一般方法不起作用。

特别是，当随机变量的值不符合分布的有效性域时，通用 nnlf 方法返回“inf”。 除非您猜测起始值非常接近最终拟合值，否则“fmin”优化器不能很好地处理此目标函数。

一般来说，.fit 方法需要对 pdf 的适用范围存在限制的分布使用约束优化器。

Answer 3

问题还在于帕累托分布的方差没有为 c 定义 < 2.

Answer 4

请参阅文档了解

fit

：

返回形状（如果适用）、位置和规模的估计值 来自数据的参数。默认估计方法是最大值 似然估计 (MLE)，但矩量法 (MM) 也是 可用

注意：

可以通过传入将一些参数固定为特定值 关键字参数

但是，您在

dataset=scipy.stats.pareto.rvs(1.5,size=10000)  #generating data

中生成了 ds，其中 shape=1.5 和默认使用的其他参数（再次参见 docs） -

 rvs(b, loc=0, scale=1, size=1, random_state=None) # Random variates

- 使用给定的

b

作为您给出的 fixed 参数它

from scipy.stats import pareto

shape= 2.7
ds = pareto.rvs(shape, size=1000)  #generating pareto distr with given ! shape
params = pareto.fit(ds, method="MLE")
print(f"shape: {params[0]:.2f}, \nmean(loc): {params[1]:.2f}, \nscale(std): {params[2]:.2f},".format(params))

##shape: 2.77,
##mean(loc): -0.00,
##scale(std): 1.00,

再次强调，- 对于任一方法（MLE 或 MM - 矩量法）：

返回的答案不保证是全局最优的；它可能 只能是局部最优，否则优化可能会完全失败

要查找参数，最好使用 curve_fit 此处使用 pareto_law_objective_func，或者可以尝试 np.polyfit