フィッティング

フィッティングの実行

実験データを解析する際に、かなりの頻度で使うことになる「フィッティング」をpytonで実行してみましょう。フィッティングとは、何らかの理論予測と実験データを照らし合わせ、実験データを最も上手く説明できるような理論予測パラメーターを推定するプロセスのことです。例を見ながら、以下で見ていきましょう。

以前使っていた例をここでも使います。 

x = np.array( [2.1, 2.9, 1.9, 3.0, 2.05, 3.1])

というnumpy array xが、実は調和振動する振り子の座標を表す実験データだとします。対応する経過時間は

t =  np.arange(x.size) * 0.1

で、表されているとしましょう。グラフにすると

のようになります。

調和振動する振り子の座標は一般に Asin⁡(2πt/T−ϕ)+BA\sin(2\pi t/T-\phi)+B と書くことができます。そこで、この理論式と実験データが一致するように、パラメーター A,T,ϕ,BA, T, \phi, B を決めることにしましょう。このとき、これらのパラメーターはフィッティングパラメーターと呼ばれます。さて、実験条件から、フィッティングパラメーターについて以下のことがわかっていたとしましょう。

  • ϕ=π/2\phi=\pi/2（値が確定している）

  • TT はおおよそ0.19

  • AA はおおよそ0.49

  • BB はおおよそ2.5

この理論予測から、まずフィッティングに必要な関数を作ります。

def pendulum(tt,period,amp,offset):
    return amp*np.sin(2*np.pi*tt/period-np.pi/2)+offset

次に、値が確定していないフィッティングパラメーターA,T,BA, T, Bの予測値を用意します。つまり

guess = [0.19, 0.49, 2.5]

のようにします（角括弧[ ]で囲まれていることに注意してください。つまりリストとして予測値をまとめています）。これで、準備は整いました。早速フィッティングしてみましょう。そのために、scipyからcurve_fit関数をインポートして、実行します。フィッティングは

from scipy.optimize import curve_fit
popt, pconv = curve_fit(pendulum, t, x, guess)

のように行います。 一行目ではscipyからcurve_fit関数をインポートしています。二行目ではcurve_fit関数に

  • フィッテングに使う関数名

  • 使う横軸データ

  • 使う縦軸データ

  • フィッテングパラメーターの推定値の入ったリスト

を渡します。するとフィッテングの結果がpoptという変数とpconvという変数に代入されます。まず、poptにはフィッテングによって求められたパラメーターの推定値がnumpy arrayとして入っています。次のセルに

popt

と実行してみれば、実際にどのようなフィッテングパラメーターの推定値が得られたかがわかります。もしも、得られた推定値が、明らかにおかしい、という場合には

  • curve_fit関数に渡すフィッティングパラメーターの予測値を変えてフィッティングをやり直す

  • そもそも考えているモデル（今回の例で言えばpendulum関数）がおかしいのではないかと考えてみる

  • 実験に不備があったかもしれないと検討してみる

といった必要があるでしょう。

次に、pconvにはフィッテングパラメーターの共分散がnumpy arrayとして入っています。詳細は省きますが、この共分散の対角成分が、フィッティングパラメーターの推定値の標準偏差を与えてくれます。より詳しくは誤差解析の教科書を読みましょう。つまり、フィッティングパラメーターの推定値の誤差の目安になります。したがって、フィッティングパラメーターの誤差を標準偏差として知りたい場合には

np.sqrt(np.diag(pcov))

を実行すれば良いことになります。

以上のフィッティングの手続きでは、実験データのそのものの誤差（測定誤差）については、気にしていませんでした。もしも、何らかの理由で実験データの誤差が分かっている場合には、それをフィッティングに反映させることができます。

もしも、各実験データについてその誤差の絶対値を標準偏差として知っている場合には、その標準偏差を格納するnumpy arrayを作り、curve_fit関数に渡します。さらにabsolute_sigmaオプションをTrueに設定します。すると、返ってくる共分散は、指定した実験データの誤差を絶対として評価して、フィッティングパラメーターの共分散を返してくれます。コードの例としては

err = np.array([0.1, 0.2, 0.15, 0.2, 0.2, 0.1])
popt, pconv = curve_fit(pendulum, t, x, guess, err, abolute_sigma=True)

のようになるでしょう。

一方で、実験データの誤差の相対値しか分かっていない場合もあると思います。例えば、1つ目の測定値の誤差が一番大きく、2つ目の測定値の誤差はその半分程度とわかっている（しかし、誤差の絶対値はわからない）というような場合です。この場合

err = np.array([[0.1, 0.2, 0.15, 0.2, 0.2, 0.1]])
popt, pconv = curve_fit(pendulum, t, x, guess, err, abolute_sigma=False)

とすることで、相対的な誤差を反映させた共分散を得ることができます。

フィッティングの可視化

せっかくフィッティングができたので、それをグラフに反映させてみましょう。まず、グラフをつくるにはnumpy arrayが必要です。フィッティング結果をうけて

fitted = pendulum(t,*popt)

を実行すると、新しいnumpy array であるfittedが作られます。ここで*はpoptというリストの中身を要素ごとにバラバラにする操作を意味しています。つまりpendulum関数に、時間を表すnumpy array tと、フィッティングパラメーターの推定値をわたしています。このデータをグラフに加えると

fig1, ax1 = plt.subplots()

curve1 = ax1.plot(t,x, marker='o',c='b')
curve2 = ax1.plot(t,fitted, marker='+', c='r')
ax1.set_xlabel("time (s)")
ax1.set_ylabel("position (m)")

のようになります。

しかし、この図では大まかにフィッティングがあっているいことはわかりますが、何も知らない人からみると、単に2つの実験結果を描いているかのようもみえます。これは

  • 実験結果のデータ点を折れ線でつないでしまっていること

  • フィッティング結果のデータ点数が少なくてなめらかでないこと

  • フィッティング結果のデータ点にマーカーがつかわれていること

に起因します。そこでこの3つをを解消しましょう。まず、フィッティング結果を元に、もっとなめらかなデータをつくります。例えば100倍なめらかなデータをつくりだければ

multifactor = 100
t_fit = np.linspace(t[0],t[-1],t.size*multifactor)

のようにします。2行目はnumpyのlinspaceという関数をつかって等差数列を作り出しています。これはarangeに似ているのですが、最初の数字と最後の数字と要素数を渡すことで等差数列を作る機能があります。つまり、と同じ範囲で100倍細かい時間データを作っていることに対応します。

今作ったデータを使ってフィッティングを可視化すると

fig1, ax1 = plt.subplots()

curve1 = ax1.plot(t,x, marker='o',c='b',linestyle='None')
curve2 = ax1.plot(t_fit,pendulum(t_fit,*popt), marker='None', c='r')
ax1.set_xlabel("time (s)")
ax1.set_ylabel("position (m)")

だいぶわかりやすくなったと思います。linestyle='None'のオプションによって、実験データから折れ線を消し、marker='None'によってフィッティング結果からマーカーを消しています。せっかくなので凡例 (legend, ラベル)も追加し、プロットの範囲も調整しましょう。

fig1, ax1 = plt.subplots()

curve1 = ax1.plot(t,x, marker='o',c='b',linestyle='None', label='exp data')
curve2 = ax1.plot(t_fit,pendulum(t_fit,*popt), marker='None', c='r', label='fitting')
ax1.set_xlabel("time (s)")
ax1.set_ylabel("position (m)")
ax1.legend(loc='upper left')
ax1.set_xlim(-0.02,0.52)
ax1.set_ylim(1.8,3.5)

最後に、もしも実験データそのものの誤差が分かっているのであれば、それもエラーバーとしてグラフに加えた方がよいかもしれません。その場合plotのかわりにerrorbar関数を使い、errorbarの大きさを指定します。

err = np.array([0.1, 0.2, 0.15, 0.2, 0.2, 0.1])
fig1, ax1 = plt.subplots()

curve1 = ax1.errorbar(t,x, yerr=err,marker='o',c='b',linestyle='None', label='exp data')
curve2 = ax1.plot(t_fit,pendulum(t_fit,*popt), marker='None', c='r', label='fitting')
ax1.set_xlabel("time (s)")
ax1.set_ylabel("position (m)")
ax1.legend(loc='upper left')
ax1.set_xlim(-0.02,0.52)
ax1.set_ylim(1.8,3.5)

最後に、フィッティングのような解析を、numpy arrayの一部にだけ適用した場合は、以前学んだスライスを使いましょう。つまり

x_part = x[0:3]

のように一部分だけを取り出したarrayを作り、フィッティングに利用すると良いでしょう。

前へグラフをもうちょっとマシにする次へデータファイルの読み書き

最終更新