測定誤差を考慮した最小二乗法

誤差を考慮した最小二乗法

理論の説明をした後に、解析的に求めた結果を反映させた関数をColabで実装します。
(大学の講義で学んだことを自分なりにかみ砕いて説明したものです)

理論的な話

普通の単回帰

状況としては、説明変数$x_i$ と、目的変数$y_i$ が実験を通して得られたという状況を考えます。
これで単回帰をすれば、

$ Y_i =ax_i+b$

という1次関数の形で、目的変数$y$ を予測するモデルが出来上がります。

そこで、

$\displaystyle \sum_i (y_i-Y_i)^2$

が最小になるような$a,b$ を見つければいいというのが普通の最小二乗法です。

誤差考慮Ver

何をするか

まず、測定をするときには誤差が生じます。その誤差が偶然誤差であり、その偶然誤差が正規分布に従う場合にのみ使用することのできるものです。導出は省きますが、結論、今回最小化したいのは

$ \chi^2=\displaystyle \sum_i \displaystyle \frac{(y_i-Y_i)^2}{\sigma_i^2} $

です。これはカイ二乗値です。
ここで、$\sigma_i$ は $y_i$ の測定誤差です。

$\chi^2$を最小化するためには、パラメータ微分が0となるような$a,b$ を探す必要があります。
これは頑張れば解析的に求めることができます。

数値的にやるのであれば勾配降下法とかを使うんだと思います。

ここでは、解析的に求める手順（細かい計算は省く）を紹介します。

微分をして、極値をとるタイミングを探す。

パラメータ微分が0となる点を見つけます。

やることとしては、

$ \displaystyle\frac{\partial \chi^2}{\partial a}=0$
$ \displaystyle\frac{\partial \chi^2}{\partial b}=0$

という連立方程式を解いて、$a,b$ を求めます。

計算すると、

$a=\displaystyle (\sum_i\frac{x_iy_i}{\sigma_i^2}\sum_i\frac{1}{\sigma_i^2}-\sum\frac{y_i}{\sigma_i^2}\sum_i\frac{x_i}{\sigma_i^2})/\Delta$

$b=\displaystyle (\sum_i\frac{x_i^2}{\sigma_i^2}\sum_i\frac{y_i}{\sigma_i^2}-\sum\frac{x_iy_i}{\sigma_i^2}\sum_i\frac{x_i}{\sigma_i^2})/\Delta$

と書くことができます。

ここで、

$\Delta=\displaystyle \sum_i\frac{x_i^2}{\sigma_i^2} \sum_i \frac{1}{\sigma_i^2}-(\sum_i\frac{x_i}{\sigma_i^2})^2$　

と置きました。

パラメータの誤差

測定誤差があれば、パラメータにも誤差があるでしょう。

これは、誤差の伝搬の考え方から

$\sigma_a=\displaystyle \sqrt{\sum_i\frac{1}{\sigma_i^2} /\Delta}$

$\sigma_b=\displaystyle \sqrt{\sum_i\frac{x_i^2}{\sigma_i^2} /\Delta}$

とすることができるそうです。

換算カイ二乗について

$ \chi_{redu}^2=\displaystyle\frac{\chi^2}{df}$

という式で換算カイ2乗を定義します。ここで、$df$は自由度(観測の長さ-パラメータ数)です。

この値が1に近いほど良い精度のモデルであるという事ができます。

実装

コードは以下の通りです。
一発で見れるように、パラメータのリストを返り値にしました。

import numpy as np

def chi2_reg(x,y,ye,Ndata,Npara=2):
  param={}  #空のパラメータリストを作成

  param['Delta'] = np.sum(x**2/ye**2) * np.sum(1/ye**2) - (np.sum(x/ye**2))**2

  param['a']=((np.sum((x*y)/ye**2))*(np.sum((1/ye**2)))-(np.sum((y/ye**2)))*(np.sum((x/ye**2))))/param['Delta']
  param['sig_a']=np.sqrt(np.sum(1/ye**2)/param['Delta'])

  param['b']=(np.sum(x**2/ye**2)*np.sum(y/ye**2)-np.sum(x*y/ye**2)*np.sum(x/ye**2))/param['Delta']

  param['sig_b']=np.sqrt(np.sum(x**2/ye**2)/param['Delta'])

  param['chi2']=np.sum(y**2/ye**2)-param['a']*np.sum(x*y/ye**2)-param['b']*np.sum(y/ye**2)
  param['red_chi2']=param['chi2']/(Ndata-Npara)

  return param

このあと、x、y、ye、Ndata=len(x)を指定し、

result=chi2_reg(x,y,ye,Ndata)
result

とすれば、結果を見ることができます。