フィッティング２(最小二乗法)

　いよいよ，準備した（疑似）測定データをつかって，最小二乗法による関数のフィッティングを行ってみよう． 　測定データには必ず誤差や雑音に由来するバラツキがあることを重々承知の上で，データ間（xとy）の関係性を見出したい． 準備したデータであれば，比例関係（直線）が見出せるはずである．

最小二乗法

以下は，（前項で準備した）測定データである（mother lineは消した）．バラツキはあるものの，いかにも比例関係がありそうである．

データ間の関係は理論的な要請などから関数形が定まっていることも多い．ここでは１次関数　\(f(x) = ax + b\)　で記述できるはずである． もしもデータにバラツキがなければ二点の情報から直線の傾きaと切片bの値を決めることができるが，残念ながらある． 言い換えると三点以上の測定データを同時に満足するa,bはない．そこで，すべての測定データとの差が最も小さくなるようにa,bを決めることにする． 最小二乗法ではその名の通り，全データについて測定データとの差の二乗の和（E）を取り，これが最小になるようにa,bを決める．

具体的には，

\[\begin{split}E &= [y_1-(ax_1+b)]^2 + [y_2-(ax_2+b)]^2 + ...... + [y_n-(ax_n+b)]^2 \\ &= \sum_{i=1}^n [y_i-(ax_i+b)]^2\end{split}\]

を最小にするaとbを探す．この場合Eは二次式なので，最小値では

\[\frac{\partial E}{\partial a} = 0, \frac{\partial E}{\partial b} = 0\]

である．偏微分をそれぞれ実行すると，以下の関係式が得られる．

\[\begin{split}\frac{\partial E}{\partial a} &= 2[y_1-(ax_1+b)](-x_1) + 2[y_2-(ax_2+b)](-x_2) + ...... + 2[y_n-(ax_n+b)](-x_n) \\ &= -\sum_{i=1}^n 2(y_i-ax_i-b)x_i = 0\end{split}\]

(1)\[\sum_{i=1}^n x_i y_i - a\sum_{i=1}^n(x_i)^2 - b\sum_{i=1}^n x_i = 0\]

\[\begin{split}\frac{\partial E}{\partial b} &= 2[y_1-(ax_1+b)](-1) + 2[y_2-(ax_2+b)](-1) + ...... + 2[y_n-(ax_n+b)](-1) \\ &= -\sum_{i=1}^n 2(y_i-ax_i-b) = 0\end{split}\]

(2)\[\sum_{i=1}^n y_i - a\sum_{i=1}^n x_i - b\sum_{i=1}^n 1 = 0\]

　(2) を (1) に代入すると，

\[a = \frac{n\sum_{i=1}^n x_i y_i - \sum_{i=1}^n x_i \sum_{i=1}^n y_i}{n\sum_{i=1}^n (x_i)^2-(\sum_{i=1}^n x_i)^2}\]

さらにこのaを (2) 式へ入れると

\[b = \frac{\sum_{i=1}^n (x_i)^2\sum_{i=1}^n y_i - \sum_{i=1}^n x_i \sum_{i=1}^n x_i y_i}{n\sum_{i=1}^n (x_i)^2-(\sum_{i=1}^n x_i)^2}\]

ということで，a,bが計算できる．

Note

（参考文献） 科学計測のためのデータ処理入門 ，河田聡編著（CQ出版，2002）

a,bの計算をpythonで行うと以下のようになる．

import numpy as np
import matplotlib.pyplot as plt

def f(x):
    return 5.*x + 10.

# Pseudo measured data by random number
xdata = np.linspace(-10, 10, num=101)
np.random.seed(1234)
ydata = f(xdata) + 5.*np.random.randn(xdata.size)

# Least squares method
x_sum = 0.
y_sum = 0.
xx_sum = 0.
yy_sum = 0.
xy_sum = 0.
for i in range(xdata.size):
    x_sum = x_sum + xdata[i]
    xx_sum = xx_sum + xdata[i]**2
    y_sum = y_sum + ydata[i]
    yy_sum = yy_sum + ydata[i]**2
    xy_sum = xy_sum + xdata[i]*ydata[i]

delta = xdata.size*xx_sum - x_sum**2
a_fit = (xdata.size*xy_sum - x_sum*y_sum)/delta
b_fit = (xx_sum*y_sum - x_sum*xy_sum)/delta
print(a_fit,b_fit)

# PLot
plt.figure(figsize=(8,5))
plt.plot(xdata,ydata,'bo', label='Exp.')
plt.plot(xdata,a_fit*xdata+b_fit,'k-', label='fitted line', linewidth=10, alpha=0.3)
plt.xlabel('x_data')
plt.ylabel('y_data')
plt.legend(loc='best',fancybox=True, shadow=True)
plt.grid(True)
plt.show()

一気にコードが長くなってしまったような気がするが，新たに書き加えたのは１２行目以降２８行目までである．

１３～１７行目：和をとるための変数の初期値を設定（かつ型を宣言）．

１８～２３行目：データの数（xdata.size）だけ和をとる．１９行目は，xdataを次々と足している，２０行目はxdataを二乗したものを足している．

２６，２７行目：aとbを算出する．

３３行目：plot内で，得られたa，bを使ってyの値を計算している．

得られたa,bの値は，

(4.9868427765835426, 10.188239290707264)

もともと測定データを作るのにつかった一次関数の傾きa=5,切片b=10であった．これにとても近い数値が得られている． グラフを描けば以下の通り．最小二乗法により得られた１次関数がデータによく沿っている．