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

冈萨雷斯，M. C.，伊达尔戈，C. A.，& Barabasi，A. L.（2008）。了解个人移动模式。 Nature，453（7196），779-782。

这样的：

the picture of equation

这是一个指数截断的幂律。有三个参数需要估计：rg0，beta和K.现在我们有几个用户的回转半径（rg），并将其上传到Github上：radius of gyrations.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来解决它。

来源

2017-04-17 Daniel Chen

请给python脚本读取数据，并指定要使用Python来解决问题。 –

看来我们需要两列变量：rg和p（rg）来估计参数，到目前为止您只给出了一列数据。你能提供更多的信息吗？ –

还请添加更多标签，例如python。 –

据@迈克尔的建议下，我们可以使用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()

来源

2017-04-18 03:41:32

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

回答

相关问题