Pythonプログラミング（ステップ７・最小二乗法・理論編・補足）

このページでは、最小二乗法で得た回帰直線のパラメータの「誤差」について考える。

理想と現実

ある量 $X$ に対する変量 $Y$ の関係が、分布が$N(0,\sigma^2)$であるような確率的な変動 $\xi$ を含んだ線形の関係 $$ Y = a X + b + \xi $$ にあるとしよう。 ただしこれは「天上」での話しであって、地上の我々は、$n$ 回のサンプリングを行って、$(x_i, y_i)$ の $n$ 組の値を得ることができるだけである。 このような状況で、最小二乗法によって推定した $a,b$ には、どれくらいのあいまいさが残されているのだろうか。

上記のモデルの「真の」パラメータ $a,b$ は $$ \begin{eqnarray} a & = & \frac{E(XY) - E(X) E(Y)}{E(XX) - E(X)^2} ,\\ b & = & \frac{E(XX) E(Y) - E(XY) E(X)}{E(XX) - E(X)^2} \end{eqnarray} $$ で与えられるが、これらの式中に現れている平均値($E(\cdot)$)が、有限のサンプルデータから得た値と、どの程度食い違うかを考えなければならない。

サンプルから得たパラメータの推定値

ここで、$y_i=a x_i + b + \xi_i$とおいて、サンプルについての平均を使って、$a$の推定値（真の値と区別するため、ここでは$\tilde{a}$と表記する）を計算すると、 $$ \begin{align} \tilde{a} & = \frac{\frac{1}{n} \sum_i x_i (a x_i + b + \xi_i) - \left( \frac{1}{n} \sum_i x_i \right) \left(\frac{1}{n} \sum (a x_i + b + \xi_i) \right)}{\frac{1}{n} \sum_i {x_i}^2 - \left(\frac{1}{n} \sum_i {x_i}\right)^2} \\ & = a + \frac{\frac{1}{n} \sum_i \left( x_i - \frac{1}{n} \sum_j x_j \right) \xi_i}{\frac{1}{n} \sum_i {x_i}^2 - \left(\frac{1}{n} \sum_i {x_i}\right)^2} \end{align} $$ となる。 $\tilde{a}$の「揺らぎ」は、$x_i$の平均値からのずれに応じた重みを付けて、揺らぎ（誤差）$\xi_i$ を線形に加算した格好になっている。

さて、正規分布$N(0,\sigma^2)$に従うような独立な確率変数 $Z_i$ を「重み」$w_i$ を付けて加えたとき、その和 $\sum_i w_i Z_i$ が従う統計性を考えよう。 それぞれの項の分散は $(w_i \sigma)^2$ になる、すなわち、 $$ w_i Z_i \sim N(0, {(w_i \sigma)}^2) $$ であるから、総和全体では $$ \sum_i w_i Z_i \sim N\left(0, \sum_i (w_i \sigma)^2\right) $$ である。 このことを使って、現実と理想の差 $\tilde{a} - a$ の分散 ${\tilde{\sigma}_a}^2$ を計算する。 上の「重み」にあたるのが $$ w_i = \frac{\frac{1}{n} \left( x_i - \frac{1}{n} \sum_j x_j \right)}{\frac{1}{n} \sum_j {x_j}^2 - \left(\frac{1}{n} \sum_j {x_j}\right)^2} $$ であることに注意すれば、直ちに、 $$ {\tilde{\sigma}_a}^2 = \frac{ \frac{1}{n} \sigma^2 }{\frac{1}{n} \sum_i {x_i}^2 - \left(\frac{1}{n} \sum_i {x_i}\right)^2} $$ が得られる。

ここで、$\xi$ の分散 $\sigma^2$ の推定値（不偏分散）は、このページの後半にあるように、 $$ \sigma^2 = \frac{1}{n-2} \sum_i \left( y_i - a x_i - b \right)^2 $$ となる。

ここまでをまとめると、データから推定されたフィッティングパラメータ $\tilde{a}$ は、「理想」の値 $a$ の周りで分散 $$ {\tilde{\sigma}_a}^2 = \frac{ 1 } { \frac{1}{n} \sum_i {x_i}^2 - \left(\frac{1}{n} \sum_i {x_i}\right)^2 } \left( \frac{\sigma^2}{n} \right) $$ だけ「ばらついて」いる可能性があり、これがパラメータ自体の曖昧さを表すひとつの指標と考えることができる。 その分母は、$X$ の分散なので、データが横軸方向に大きく散布していれば、曖昧さは減じられる。 一方、縦軸方向に回帰直線から大きく分散していれば、当然ながら、曖昧さは増加する。 また、分散が同じでも、データ数 $n$ が大きいほど、曖昧さは小さくなる。

全く同様に考えることで、y切片 $\tilde{b}$についても、その分散 ${\tilde{\sigma}_b}^2$ を見積もることができて、結果は $$ {\tilde{\sigma}_b}^2 = \frac{\frac{1}{n} \sum_i {x_i}^2 } { \frac{1}{n} \sum_i {x_i}^2 - \left(\frac{1}{n} \sum_i {x_i}\right)^2 } \left( \frac{\sigma^2}{n} \right) $$ となる。 こちらは、$X$ の平均が0に近いほど分散が小さくなる、という形の式になっている。

分散$\sigma^2$の推定についての考察

この節では、$\xi$の母分散の推定値について考察する。 我々は、元々 $a,b,\sigma^2$ についての知識を持たず、データ$\{(x_i,y_i)\}$からそれを推定しようとしているわけだ。 これを、ベイズの定理を使って表すと $$ P\left(a,b,\sigma^2 \middle| \{(x_i,y_i)\} \right) = \frac{P\left( \{(x_i,y_i)\} \middle| a,b,\sigma^2 \right) P(a,b,\sigma^2)}{P\left(\{(x_i,y_i)\}\right)} $$ ということになる。

ここで、右辺の分母はデータが与えられれば「定数」と見なせるので、$a, b, \sigma^2$の推定には直接関係しない。 一方、右辺の分子の「事前確率」$P(a,b,\sigma^2)$ は、我々が$a,b,\sigma^2$の値の組をどれくらい「ありそう」に思っているかの程度（degree of belief)を表している。 どんな値でも「あり」であれば $P(a,b,\sigma^2)=\textrm{const.}$ と考えるのが自然であろう。 ところが、「モデルの動作域」まで考えると、事情はそれほど単純ではない。

データの確率分布を $a,b,\sigma^2$ の関数として眺めると、それは $$ P\left( \{(x_i,y_i)\} \middle| a,b,\sigma^2 \right) = \left( \frac{1}{\sqrt{2\pi \sigma^2}} \right)^n \exp\left( -\frac{F(a,b)}{2\sigma^2} \right) $$ のような格好をしていた。 $a,b$ 平面で、$F(a,b)$ が最小となる点を基準に等高線を描くと、最小二乗法で決まる最小値の周りを囲むような（主軸が傾いた）楕円形になっている。 ここで、$\sigma$ が大きくなると等高線を示す楕円は相似的に（長さが $\sigma$ 倍に）大きくなり、小さくなると楕円も縮む。 ただし、最小値を与える $(a,b)$ の位置は変わらない。

このことを視点を変えて見ると、パラメータ $a,b$ の「実質的に動ける面積」が $\sigma^2$ に比例して変化する、とも言える。

ベイズの定理に立ち戻ると、事前確率に対して $$ P(a,b,\sigma^2) = C \sigma^2 $$ ($C$は定数)のような重み付けをしておかないと、パラメータ $a,b$ が等しい重みで「ありそう」と考えていることにならない （さもないと、$\sigma$によって$(a,b,\sigma^2)$の各点での事前確率の大きさが変わる）、ということになる。

事前確率まで考慮して、尤度 $P\left(a,b,\sigma^2 \middle| \{(x_i,y_i)\} \right)$ を最大化するような $\sigma^2$ を求めると、 $$ \sigma^2 = \frac{1}{n-2} \sum_i^n \left(y_i - a x_i - b \right)^2 $$ となる。

上記の議論をさらに一般化し、パラメータの次元を増やして、超平面 $y = a_1 x_1 + a_2 x_2 + \cdots + a_{q-1} x_{q-1} + b$ で最小二乗法を行ったとしよう（この場合、パラメータ空間の次元は $q$）。 すると、 $P(a_1,\cdots,a_{q-1},b,\sigma^2) \sim \sigma^q$ となり、 上式の $n-2$ のところは、$n-q$ となる。

上記は、統計学の教科書などに記されている、統計的な自由度まで考慮した標準誤差（standard error of estimate)に一致する。