如何使用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估計指數截斷冪律的參數？

回答

相關問題