物理の駅 Physics station by 現役研究者

テクノロジーは共有されてこそ栄える

Pythonでscipyを使わず一次関数(y=ax+b)でフィッティングし係数と係数の誤差を得る

参考にした解説:

http://www.cc.u-ryukyu.ac.jp/~fukami/p0.pdf

及び岐阜大学教育学部の物理学実験及びコンピュータ処理の解答例

ある二次元(x,y)の点列で与えられるデータ群を、一次関数 (=直線)

 \displaystyle
y=a x + b

でフィッティングするとき、カイ二乗を最小にすればよい。yにつく誤差が同じだとすると最小二乗法で解けて、残差

 \displaystyle
S = \sum(ax_i+b-y_i)^{2}

を最小にするのと同等になる。それぞれの係数の偏微分を\frac{\partial S}{\partial a}, \frac{\partial S}{\partial b}を0にすればよいので、

 \displaystyle
a = \frac{N\sum x_i y_i - \sum x_i \sum y_i}{N\sum x^{2}_i - (\sum x_i)^{2}},

 \displaystyle
b = \frac{\sum x^{2}_i \sum y_i - \sum x_i \sum x_i y_i}{N\sum x^{2}_i - (\sum x_i)^{2}}

と得ることが出来る。

係数2つの自由度を加味したyの誤差\sigma_y

 \displaystyle
\sigma_y = \sqrt{\frac{1}{N-2}\sum(ax_i+b-y_i)^{2}}

となる。

誤差伝搬により、係数 a, bの誤差は、

 \displaystyle
\sigma_a = \sigma_y \sqrt{\frac{N}{N\sum x^{2}_i-(\sum x_i)^{2}}}

 \displaystyle
\sigma_b = \sigma_y \sqrt{\frac{\sum x^{2}_i}{N\sum x^{2}_i - (\sum x_i)^{2}}}

となる。

さて、これをPythonで実装してみる。

import numpy as np

# 再現性を確保するために乱数のSEEDを固定
np.random.seed(1)

# データ生成
x = np.arange(6, dtype=float)
y = 3*x + 2
# 乱数で誤差を付与する
x += np.random.normal(loc=0, scale=0.1, size=6)
y += np.random.normal(loc=0, scale=0.1, size=6)

# 頻出する総和を先に計算
N   = len(y)
Nxy = np.sum([xi*yi for xi,yi in zip(x,y)])
Nx  = np.sum([xi    for xi,yi in zip(x,y)])
Ny  = np.sum([yi    for xi,yi in zip(x,y)])
Nx2 = np.sum([xi*xi for xi,yi in zip(x,y)])

# 係数
a = (N   * Nxy - Nx * Ny )/(N * Nx2 - Nx**2)
b = (Nx2 * Ny  - Nx * Nxy)/(N * Nx2 - Nx**2)

# Yの誤差
sigma_y = np.sqrt(1/(N-2) * np.sum([(a*xi+b-yi)**2 for xi,yi in zip(x,y)]))

# 係数の誤差
sigma_a = sigma_y * np.sqrt(N   / (N * Nx2 - Nx**2))
sigma_b = sigma_y * np.sqrt(Nx2 / (N * Nx2 - Nx**2))

# グラフ生成
import matplotlib.pyplot as plt

xx=np.linspace(0,5,100)
yy=a*xx+b

plt.plot(x,y,'ro',label='data')
plt.plot(xx,yy,label="y = {:.3f} $\pm$ {:.3f} x + {:.3f} $\pm$ {:.3f}\nsigma$_y$={:.3f}".
         format(a,sigma_a,b,sigma_b,sigma_y))
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

上記の数式を愚直に実装したが、numpyにはN次元関数でフィッティングする関数が用意されているので、もっと簡単に記述できる。 係数の誤差は、分散共分散行列 {\rm cov} で与えられる。


{\rm cov}=
\begin{pmatrix}
\sigma_a^{2} & \sigma_a \sigma_b \\\
\sigma_b \sigma_a & \sigma_b^{2}
\end{pmatrix}

この、対角成分だけ使って、係数の誤差 \sigma_a, \sigma_bを計算した。

import numpy as np

# 再現性を確保するために乱数のSEEDを固定
np.random.seed(1)

# データ生成
x = np.arange(6, dtype=float)
y = 3*x + 2
# 乱数で誤差を付与する
x += np.random.normal(loc=0, scale=0.1, size=6)
y += np.random.normal(loc=0, scale=0.1, size=6)
N   = len(y)

# np.polyfitで係数と係数の誤差を計算
p, cov =np.polyfit(x,y,1, cov=True)
a = p[0]
b = p[1]
sigma_a = np.sqrt(cov[0,0])
sigma_b = np.sqrt(cov[1,1])

# Yの誤差
sigma_y = np.sqrt(1/(N-2) * np.sum([(a*xi+b-yi)**2 for xi,yi in zip(x,y)]))

# グラフ生成
import matplotlib.pyplot as plt

xx=np.linspace(0,5,100)
yy=np.polyval(p,xx) #結果からyを計算

plt.plot(x,y,'ro',label='data')
plt.plot(xx,yy,label="y = {:.3f} $\pm$ {:.3f} x + {:.3f} $\pm$ {:.3f}\nsigma$_y$={:.3f}".
         format(a,sigma_a,b,sigma_b,sigma_y))
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()

3次関数でも4次関数でも、全て線形なので同様に解ける。

誤差付きのデータの場合は以下のようにすればよい

import numpy as np

# データ生成
x = np.array([0,1,2,3,4,5])
y = np.array([0,4,6,8,10,15])
yerr = np.array([0.4,1,1,1,1,0.4])
N   = len(y)

# np.polyfitで係数と係数の誤差を計算.誤差はw=1/sigmaで与える
p, cov =np.polyfit(x,y,1, w=1.0/yerr,cov=True)
a = p[0]
b = p[1]
sigma_a = np.sqrt(cov[0,0])
sigma_b = np.sqrt(cov[1,1])

# グラフ生成
import matplotlib.pyplot as plt

xx=np.linspace(0,5,100)
yy=np.polyval(p,xx) #結果からyを計算

plt.errorbar(x,y,yerr=yerr,fmt='o',label='data',capsize=5,ecolor='black',markeredgecolor="black",color='w')
plt.plot(xx,yy,label="y = ({:.3f} $\pm$ {:.3f}) x + ({:.3f} $\pm$ {:.3f})".format(a,sigma_a,b,sigma_b))
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.grid()
plt.show()

メモ: 線形回帰 line fit