Linear Regression and Uncertainty

観測した \(N\) 個のデータの組に \(\hat{y}_i = ax_i +b\) という線形の関係があるとします. パラメタ \(a\), \(b\) を最小二乗法によって決定することを考えます. ここで観測誤差 \(e_i\) がパラメタ \(a\), \(b\) の精度に与える影響を定量化します. パラメタのばらつき \(\sigma_a\), \(\sigma_b\) を総データ数 \(N\), 説明変数のばらつき \(\sigma_x\), 説明変数の平均値 \(\bar{x}\), 観測誤差のばらつき \(\sigma_e\) で表現するための手続きを示します.

線形回帰

1 次線形回帰の公式

\(N\) 個の観測量ペア \((x_i,y_i)\) に対して \(\hat{y}_i = ax_i +b\) という線形モデルでフィッティングを実行します. 最小二乗法によって \(a\), \(b\) は以下の式で与えられます.

\[ a = \frac{N\sum_i y_ix_i - \sum_iy_i\sum_ix_i} {N\sum_ix_i^2 - \left(\sum_i x_i\right)^2}, \quad b = \frac{\sum_iy_i\sum_ix_i^2 - \sum_ix_i\sum_iy_ix_i} {N\sum_ix_i^2 - \left(\sum_i x_i\right)^2}, \]

パラメタの誤差推定

ここで, 実際の観測値には \(y_i = ax_i + b + e_i\) と誤差がのっているとします. 誤差 \(e_i\) は \(E(e) = 0\), \(V(e) = \sigma_e^2\) を満たすような相関のないノイズであるとします. このノイズ成分の大きさ \(\sigma_e\) がパラメタ \(a\), \(b\) を決定する精度に与える影響を調べます.

説明変数 \(x\) の平均値を \(\bar{x}\, (\equiv E(x))\) とします. 簡単のために \(x \leftarrow x - \bar{x}\) として \(\bar{x} = 0\) となるように変数変換します.

\[ a = \frac{\sum_i y_ix_i}{\sum_ix_i^2}, \quad b' = \frac{\sum_iy_i\sum_ix_i^2}{N\sum_ix_i^2}, \]

パラメタ \(a\) の分散を計算すると以下のようになります.¹

\[ \sigma_a^2 = \frac{\sum_i V\left(y_ix_i\right)}{\left(\sum_i x_i^2 \right)^2} = \frac{N\sigma_e^2}{(N-1)^2\sigma_x^2}, \]

ただし \(\sigma_x\) は \(x\) の分散です. パラメタ \(b'\) の分散は以下のようになります.

\[ \sigma_b^2 = \frac{V\left(\sum_i y_i \sum_i x_i^2 \right)}{\left( N\sum_i x_i^2 \right)^2} = \frac{N\sigma_e^2}{(N-1)^2}, \]

本来の切片 \(b\) は \(b = b' - a\bar{x}\) という関係にあります. \(\bar{x}\) は定数です. ここで \(b'\) と \(a\) の分散には相関がないとして² \(b\) の分散を計算すると以下の式を得ます.

\[ \sigma_b^2 = \frac{N\sigma_e^2}{(N-1)^2} + \sigma_a^2\bar{x}^2, \]

よってパラメタ \(a\), \(b\) の標準偏差はそれぞれ以下のように導出できます.

\[ \sigma_a = \sqrt{\frac{N\sigma_e^2}{(N-1)^2\sigma_x^2}}, \quad \sigma_b = \sqrt{\frac{N\sigma_e^2}{(N-1)^2} + \sigma_a^2\bar{x}^2}. \]

比例係数の誤差は \(\sigma_a \sim \sigma_e/\sigma_x\) の関数になっています. データ数 \(N\) を増やすとおよそ \(\sqrt{N}\) の関係で精度が向上することも分かります. 切片の誤差は \(\bar{x} = 0\) であれば \(\sigma_b \sim \sigma_e/\sqrt{N}\) の精度で決まる. \(\bar{x} \not= 0\) の場合には, 傾きの不定性が切片の不定性に影響を与えます. \(\bar{x}\) に比例するような項が誤差として効いてきます.

シミュレーションによる確認

\(a = 1.1\), \(b = 100\), \(\sigma_e = 21\), \(x \sim {\mathsf U}(0,100)\) (一様分布, \(\sigma_x \sim 28.9\), \(\bar{x} = 50\)) としてモンテカルロ法によって \(10^4\) 回の試行で推定したパラメタ \(a\), \(b\) のヒストグラムを作成しました. 図 1 の青い階段がシミュレーションによって求めたヒストグラムです. 解析的に推定したパラメタのばらつきを \(\sigma_a\), \(\sigma_b\) に対応する Gaussian 分布として緑色の実線で示しました. 解析的な推定によってシミュレーション結果がよく再現されています.

\(\bar{x}=0\) となるように変換することは, データ全体を \(x\) 軸方向にシフトするという操作に対応します. どちらの座標系で計算しても \(a\) についての不定性は同じ値になります. ↩
\(\bar{x}=0\) の場合には回帰直線は \(a\) の値に依らず \((0,b')\) という点を通ります. そのため \(a\) の不定性と \(b'\) の不定性は独立だとみなしても良さそうです. ↩

Last update: April 22, 2024