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\) は以下の式で与えられます.
パラメタの誤差推定
ここで, 実際の観測値には \(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\) の分散を計算すると以下のようになります.1
ただし \(\sigma_x\) は \(x\) の分散です. パラメタ \(b'\) の分散は以下のようになります.
本来の切片 \(b\) は \(b = b' - a\bar{x}\) という関係にあります. \(\bar{x}\) は定数です. ここで \(b'\) と \(a\) の分散には相関がないとして2 \(b\) の分散を計算すると以下の式を得ます.
よってパラメタ \(a\), \(b\) の標準偏差はそれぞれ以下のように導出できます.
比例係数の誤差は \(\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 分布として緑色の実線で示しました. 解析的な推定によってシミュレーション結果がよく再現されています.