Pythonを使用して指数関数の法則のパラメータを推定する方法は？以下の記事は指数関数的に切り捨て指数法則の方程式があり

Question

：

ゴンザレス、M. C.、イダルゴ、C. A.、& Barabasi、A. L.（2008）。個々の人間のモビリティパターンを理解する。 Nature、453（7196）、779-782。このような

：

それが指数切り捨てべき乗則です。 rg0、beta、Kの3つのパラメータがあります。ユーザーの回転半径（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() 

は、どのように私は上記の3つのパラメータを推定するために、RGSのこれらのデータを使用することができますか？私はPythonを使用してそれを解決することを望む。

Answer 1

@Michaelさんの提案によると、我々は問題を解決することができた結果を以下に示す。

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 

scipy.optimize.curve_fitを使用して：

[ 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()