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