物理の駅 by onsanai

Physics Station → PhSt 質問・疑問・間違いの指摘は、コメントに書くか、直接伝えるときっと良いことがあります。主にWindows or Ubuntu用の記事です

Python: 二次元データを一次関数(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()

f:id:onsanai:20200123111654p:plain

上記の数式を愚直に実装したが、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)

# 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=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()

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

メモ: 線形回帰