Python中的幂律分布拟合答案

【问题标题】：Power law distribution fitting in PythonPython中的幂律分布拟合
【发布时间】：2020-04-02 15:20:07
【问题描述】：

我正在使用不同的 python 来拟合数据集的密度函数。该数据集由从 1 秒开始的正时间值组成。

我测试了来自scipy.statistics 和powerlaw 库的不同密度函数，以及使用scipy.optimize 的函数curve_fit() 的我自己的函数。

到目前为止，我在拟合以下“修改”幂律函数时获得了最好的结果：

def funct(x, alpha, x0):
    return((x+x0)**(-alpha))

我的代码如下：

bins = range(1,int(s_distrib.max())+2,1)
y_data, x_data = np.histogram(s_distrib, bins=bins, density=True)
x_data = x_data[:-1]

param_bounds=([0,-np.inf],[np.inf,np.inf])
fit = opt.curve_fit(funct,
                    x_data,
                    y_data,
                    bounds=param_bounds) # you can pass guess for the parameters/errors
alpha,x0 = fit[0]
print(fit[0])

C = 1/integrate.quad(lambda t: funct(t,alpha,x0),1,np.inf)[0]

# Calculate fitted PDF and error with fit in distribution
pdf = [C*funct(x,alpha,x0) for x in x_data]
sse = np.sum(np.power(y_data - pdf, 2.0))
print(sse)

fig, ax = plt.subplots(figsize=(6,4))
ax.loglog(x_data, y_data, basex=10, basey=10,linestyle='None',  marker='.')
ax.loglog(x_data, pdf, basex=10, basey=10,linestyle='None',  marker='.')

拟合返回 x0 的值 8.48 和 alpha 的 1.40。在 loglog 图中，数据和拟合图如下所示：

我的第一个问题是技术问题。为什么在 funct 函数中将 (x+x0) 更改为 (x-x0) 时，opt.curve_fit 会出现以下警告和错误？由于我对 x0 的限制是 (-inf, +inf)，我希望拟合返回 -8.48。

/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:3: RuntimeWarning: 除以零遇到倒数这与 ipykernel 包是分开的，所以我们可以避免导入，直到 ValueError：残差在初始点不是有限的。

我的其他问题是理论上的。 (x+x0)^(-alpha) 是标准分布吗？ x0值代表什么，这个8.48s值如何物理解释？据我了解，这意味着我的分布对应于转移的幂律分布？我可以认为 x0 对应于将数据拟合到幂律时通常需要的 xmin 值吗？
关于这个 xmin 值，我知道在拟合过程中只考虑大于这个阈值的数据来表征分布的尾部是有意义的。但是，我想知道用 xmin 之后的幂律分布和 xmin 之前的其他分布来表征完整数据的标准方法是什么。

这是很多问题，因为我对这个主题非常不熟悉，任何评论和回答，即使是部分的，将不胜感激！

【问题讨论】：

标签： python distribution power-law scipy.stats

【解决方案1】：

(x+x0)^(-alpha) 是标准分布吗？

回答你的第二个问题，是的，它是标准发行版，称为Zipf distribution。它在 Python/NumPy as well 中实现。

x0值代表什么

这是移位参数。标准参数（如 Zipf 中的功率参数）之上的任何分布都可能具有移位和比例参数，这基本上表示您的 X 值是用不同的单位和不同的原点测量的。

关于这个 xmin 值，我知道在拟合过程中只考虑大于这个阈值的数据来表征分布的尾部是有意义的。

这是 Zipf 定律的定义方式，从 0 到无穷大。改变它意味着你的出身会不同

【讨论】：