参考にした解説:
http://www.cc.u-ryukyu.ac.jp/~fukami/p0.pdf
及び岐阜大学教育学部の物理学実験及びコンピュータ処理の解答例
ある二次元(x,y)の点列で与えられるデータ群を、一次関数 (=直線)
でフィッティングするとき、カイ二乗を最小にすればよい。yにつく誤差が同じだとすると最小二乗法で解けて、残差
を最小にするのと同等になる。それぞれの係数の偏微分をを0にすればよいので、
と得ることが出来る。
係数2つの自由度を加味したの誤差は
となる。
誤差伝搬により、係数 の誤差は、
となる。
さて、これをPythonで実装してみる。
import numpy as np
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)
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次元関数でフィッティングする関数が用意されているので、もっと簡単に記述できる。
係数の誤差は、分散共分散行列 で与えられる。
この、対角成分だけ使って、係数の誤差 を計算した。
import numpy as np
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)
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])
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)
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)
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)
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