2017-12-11 42 views
1

Python lmfitを使用して、2005-2016年の月間平均データで最小自乗近似を行います。フーリエ項が機能しなかったようだPython lmfitでカーブフィッティングするためのパラメータ推定

# t is in fractional years, e.g. 2017+122./365. 
def fun(t, a, b, c, A1, A2, A3, A4, B1, B2, B3, B4): 
    An=[A1,A2,A3,A4] 
    Bn=[B1,B2,B3,B4] 
    sum=np.sum([An[i] * np.sin(2 * np.pi * (i + 1) * t+Bn[i]) for i in range(len(An))]) 
    return a+b*t+c*t*t+sum 

mod = Model(fun) 
pars = mod.make_params(a=-10, b=0.003, c=0.01, A1=-1., A2=1., A3=1., A4=1., B1=-1., B2=1., B3=1., B4=1.) 

result = mod.fit(y, pars, t=t) 
print(result.fit_report()) 

plt.plot(t, y, 'bo') 
plt.plot(t, result.best_fit, 'r-') 
plt.show() 

fitted line and the original data dotsequation と下図のように元のコード:私は以下のような機能を構築しました。したがって、私は、A1,A2,A3 ...のような関数パラメータの適切な初期推定をどのように与えるかが不思議です。

+3

ノートであなたの

sum=np.sum([An[i]*np.sin(2*np.pi*(i+1)*t+Bn[i]) for i in range(len(An))]) 

を交換してみてください。 – mikuszefski

+1

さらに、最小限の作業例、つまりすべての 'import'を含む少なくとも完全な作業コードを提供する必要があります。最小限のデータセットですばらしいこともあります。 – mikuszefski

+0

それとは別に、 'a'は' mean'から確認できます。 'B 'はおそらく最高にゼロに設定されています。 'A1 *は0.5 *(max() - min())のようなものになります。他の「A」は0になります。 (すべてのデータがリンクに表示されているものと似ていると仮定します)、 'b'と' c'もゼロにします。 – mikuszefski

答えて

2

np.sumあなたがしたいことはしません。式は、tと同じ長さの配列ではなく、単一のスカラー値に合計されます。そのスカラー値は、次にパラメータA1、... B4を単一の値に折りたたみ、これらの値を決定する方法はありません。

形状の2次元配列(4、len(t))を作成し、最初の次元だけを合計し、len(t)の配列を4フーリエ成分に渡って残したいと思います。変数名として `sum`を使用しない良いアイデア:

sum=np.array([An[i]*np.sin(2*np.pi*(i+1)*t+Bn[i]) for i in range(len(An))]).sum(axis=0) 
+0

にあります。あなたの役に立つ答えをありがとう!実際には、私は合計条件の形に気付かなかった。 あなたの提案通りに機能しました。 https://groups.google.com/forum/#!topic/lmfit-py/u413sR1XYV0 – Cheng

関連する問題