python - 如何使用 Python 估计指数截断幂律的参数？

Question

下面的文章中有一个指数截断幂律方程：

Gonzalez, MC, Hidalgo, CA, & Barabasi, AL (2008)。了解个人的人类流动模式。自然，453（7196），779-782。

像这样：

它是一个指数截断幂律。需要估计三个参数：rg0、beta和K。现在我们得到了几个用户的回转半径（rg），并上传到Github：回转半径.txt

以下代码可用于读取数据并计算 P(rg)：

import numpy as np
# read radius of gyration from file
rg = []
with open('/path-to-the-data/radius of gyrations.txt', 'r') as f:
    for i in f:
        rg.append(float(i.strip('\n')))
# calculate P(rg)
rg = sorted(rg, reverse=True)
rg = np.array(rg)
prg = np.arange(len(sorted_data)) / float(len(sorted_data)-1)

或者您可以直接获取 rg 和 prg 数据，如下所示：

rg = np.array([ 20.7863444 ,   9.40547933,   8.70934714,   8.62690145,
     7.16978087,   7.02575052,   6.45280959,   6.44755478,
     5.16630287,   5.16092884,   5.15618737,   5.05610068,
     4.87023561,   4.66753197,   4.41807645,   4.2635671 ,
     3.54454372,   2.7087178 ,   2.39016885,   1.9483156 ,
     1.78393238,   1.75432688,   1.12789787,   1.02098332,
     0.92653501,   0.32586582,   0.1514813 ,   0.09722761,
     0.        ,   0.        ])

 prg = np.array([ 0.        ,  0.03448276,  0.06896552,  0.10344828,  0.13793103,
    0.17241379,  0.20689655,  0.24137931,  0.27586207,  0.31034483,
    0.34482759,  0.37931034,  0.4137931 ,  0.44827586,  0.48275862,
    0.51724138,  0.55172414,  0.5862069 ,  0.62068966,  0.65517241,
    0.68965517,  0.72413793,  0.75862069,  0.79310345,  0.82758621,
    0.86206897,  0.89655172,  0.93103448,  0.96551724,  1.        ])

我可以使用以下 python 脚本绘制 P(r_g) 和 r_g：

import matplotlib.pyplot as plt
%matplotlib inline

plt.plot(rg, prg, 'bs', alpha = 0.3)
# roughly estimated params:
# rg0=1.8, beta=0.15, K=5
plt.plot(rg, (rg+1.8)**-.15*np.exp(-rg/5))
plt.yscale('log')
plt.xscale('log')
plt.xlabel('$r_g$', fontsize = 20)
plt.ylabel('$P(r_g)$', fontsize = 20)
plt.show()

如何使用 rgs 的这些数据来估计上面的三个参数？我希望使用python解决它。

Answer 1

根据@Michael 的建议，我们可以使用scipy.optimize.curve_fit解决问题

def func(rg, rg0, beta, K):
    return (rg + rg0) ** (-beta) * np.exp(-rg / K)

from scipy import optimize
popt, pcov = optimize.curve_fit(func, rg, prg, p0=[1.8, 0.15, 5])
print popt
print pcov

结果如下：

[  1.04303608e+03   3.02058550e-03   4.85784945e+00]

[[  1.38243336e+18  -6.14278286e+11  -1.14784675e+11]
 [ -6.14278286e+11   2.72951900e+05   5.10040746e+04]
 [ -1.14784675e+11   5.10040746e+04   9.53072925e+03]]

然后我们可以通过绘制拟合曲线来检查结果。

%matplotlib inline
import matplotlib.pyplot as plt

plt.plot(rg, prg, 'bs', alpha = 0.3)
plt.plot(rg, (rg+popt[0])**-(popt[1])*np.exp(-rg/popt[2]) )
plt.yscale('log')
plt.xscale('log')
plt.xlabel('$r_g$', fontsize = 20)
plt.ylabel('$P(r_g)$', fontsize = 20)
plt.show()