2015-11-17 59 views
6

我試圖將一些數據擬合成指數截斷的冪律函數。我用numpy生成一些數據,我試圖用scipy.optimization來適應這些數據。 這裏是我的代碼:使用scipy的錯誤擬合curve_fit

import numpy as np 
from scipy.optimize import curve_fit 

def func(x, A, B, alpha): 
    return A * x**alpha * np.exp(B * x) 

xdata = np.linspace(1, 10**8, 1000) 
ydata = func(xdata, 0.004, -2*10**-8, -0.75) 
popt, pcov = curve_fit(func, xdata, ydata) 
print popt 

我得到的結果是:1,1,1],它不與數據相對應。 ¿我做錯了什麼?

+1

我的方法休息,因爲你的最終數據點(在655642210)快圍繞20倍值爲0.當你記錄它時,你會得到一個NaN。我已經使用我的方法計算了適合度,排除這一點並得到看起來合理的東西A = 0.00326,alpha = -0.767,B = -1.88e-8 –

+0

是的你是對的!我注意到了我的錯誤,並刪除了第二個提示。非常感謝 – ivangtorre

回答

4

雖然XNX給你答案,爲什麼curve_fit失敗在這裏,我認爲我會d提出了一種解決不依賴梯度下降(因此是合理的初始猜測)的功能形式擬合問題的不同方式。

請注意,如果你把日誌,你是裝修你的形式

\log f = \log A + \alpha \log x + B x

這是線性的每一個未知參數的功能(日誌A,α,B)

我們因此可以使用線性代數的機械解決這個由如在矩陣的形式寫方程

日誌Y =量m p

其中log y是記錄您YDATA點中的一個列向量中,p是未知參數的列向量,M爲[[1], [log x], [x]]

或明確

enter image description here

最佳擬合參數矢量然後可有效地通過使用np.linalg.lstsq

你的例子中發現的基質代碼中的問題可以寫成:

import numpy as np 

def func(x, A, B, alpha): 
    return A * x**alpha * np.exp(B * x) 

A_true = 0.004 
alpha_true = -0.75 
B_true = -2*10**-8 

xdata = np.linspace(1, 10**8, 1000) 
ydata = func(xdata, A_true, B_true, alpha_true) 

M = np.vstack([np.ones(len(xdata)), np.log(xdata), xdata]).T 

logA, alpha, B = np.linalg.lstsq(M, np.log(ydata))[0] 

print "A =", np.exp(logA) 
print "alpha =", alpha 
print "B =", B 

,恢復很好的初始參數:

A = 0.00400000003736 
alpha = -0.750000000928 
B = -1.9999999934e-08 

另外請注意,這種方法比使用curve_fit爲手頭的問題

In [8]: %timeit np.linalg.lstsq(np.vstack([np.ones(len(xdata)), np.log(xdata), xdata]).T, np.log(ydata)) 
10000 loops, best of 3: 169 µs per loop 


In [2]: %timeit curve_fit(func, xdata, ydata, [0.01, -5e-7, -0.4]) 
100 loops, best of 3: 4.44 ms per loop 
2

顯然,您的初始猜測(默認爲[1,1,1],因爲您沒有給出 - 請參閱the docs)距離實際參數太遠以使算法收斂。主要的問題可能是B,如果是肯定的,它會將您的指數函數發送到您提供的xdata的非常大的值。

嘗試提供更接近於實戰的參數一點,它的工作原理:

p0 = 0.01, -5e-7, -0.4 # Initial guess for the parameters 
popt, pcov = curve_fit(func, xdata, ydata, p0) 
print popt 

輸出:

[ 4.00000000e-03 -2.00000000e-08 -7.50000000e-01]