2013-06-06 29 views
0

我試圖使用Leastsq來擬合一條非常簡單的曲線。但是,其解決方案並未得到優化。任何人都可以給我一些建議嗎? 下面是我的代碼:SciPy LeastSq無法提供最佳解決方案

from scipy import optimize 
import numpy as np 

hl_obs = np.array([10.0, 23.0, 20.0]) 
ph=np.array([5.0,7.0,9.0]) 

tp=60 

def residuals(k_abn, ph, tp, hl_obs): 
     hr=np.log(2)/(hl_obs*24.0) 
     ph_adj=6013.79/(tp+273.15) + 23.6521*np.log10(tp+273.15)-64.7013 
     err = peval(k_abn, ph, ph_adj)-hr 
     return err 

def peval(k_abn, ph, ph_adj): 
    temp= k_abn[0]*np.power(10,-ph) + k_abn[1] + k_abn[2]*np.power(10,(-ph_adj + ph)) 
    return temp 

k_abn =np.array([1, 0, 0]) 

from scipy.optimize import leastsq 
p,ier = leastsq(residuals, k_abn, args=(hl_obs, ph, tp), maxfev=2000000) 
print p, ier 

從Excel的求解器,我知道解決的辦法應該是k_abn=[165, 0.0, 2.14]。有一次,我喂Excel的解決function peval,它產生正確的答案...

peval([165,0.0,2.14], 5.0, 13.01573)=0.002888113 
peval([165,0.0,2.14], 7.0, 13.01573)=0.001255701 
peval([165,0.0,2.14], 9.0, 13.01573)=0.001444057 

此外,我試圖增加使用epsfcn=np.finfo(np.float32).eps 任何人都可以給我一些建議精度?謝謝!

回答

2

如果你按照正確的順序ARGS:

p,ier = leastsq(residuals, k_abn, args=(ph, tp, hl_obs), maxfev=2000000) 

你得到的結果:

[ 1.65096852e+02 1.23712405e-03 2.14392540e+00] 
+0

感謝您的幫助。我很蠢...... –

+1

不傻。這是進展情況。至少在我自己的代碼中,修復1 + 1 = 3個bug總是最重要的。一旦我完成了,更復雜的(通常)自己照顧自己。它*值得花一點時間思考你可以開發哪些實踐來防止這樣的錯誤在未來發生。在這種情況下它不起作用(因爲leastsq不會讓你),但傳遞命名參數ph = ph,tp = tp等有時甚至是多餘的,因爲它可以防止這樣的錯誤發生。 – Rick