最小二乗法と最尤推定

1年近く下書きに眠っていたのを公開。これは統計的手法を真面目に勉強しはじめて、最初に「面白いな」と思った事柄だから、思い出深い。最小二乗法という昔から知っていた手法が、最尤推定という自分にとって目新しい知識とつながる体験は快かった。

回帰とは、モデルからの予測値と実測値のズレ(誤差, 残差)を最小化するようにフィッティングすることである。

具体的にいうと、説明変数(独立変数ともいう)の列x_iがあって、目的変数の列y_iが対応している。我々はモデル f(x) を用意して、モデルからの予測値、すなわちモデルに説明変数を突っ込んだ時の出力 f_i = f(x_i) がなるべくy_iに近くなるようにしたいのである。
モデルといっても、その範囲は広い。極論すれば、x_i に対して y_i を返すような関数も立派なモデルである。しかしそれは実用上何の訳にも立たないので、多くの場合は、より少ないパラメータ \mathbf{\theta} によって調節可能なモデル f(x, \mathbf{\theta})を考えるのである。いくつのパラメータを使うべきか、どのようなモデル化が望ましいかは、モデル選択といって奥が深い問題だ。いずれ考察しよう。

さて、モデルからの予測値 f_i と実測値とのズレをどう測るかには、いろいろな方法が可能である。例えば、差の絶対値の和 \sum|f_i - y_i|を取るのもありだろう。実際によく使われているのは、残差の二乗和 \sum(f_i - y_i)^2 を使う方法であり、この場合の回帰を最小二乗法と呼んでいる。

ここまでは漠然と y_i \approx f_i と考えていた。ここで発想を変えよう。モデルが系の真のあるべき姿を表していて、予測値のほうが「完璧」だと考えるのである。そこから実測値がズレてしまうのは、測定機器やら何やらの誤差であると。この場合は、y_i = f_i + \epsilon_i という具合に近似記号が等号に変わる。誤差(ノイズ) は \epsilon_i で表現されている。変形して \epsilon_i = y_i - f_i とする。誤差は確率変数であり、なんらかの分布に従うと考えるのが新しい発想法だ。

誤差というからには、平均が 0 で 0 に最大値を持つような分布であるのが自然である。ここで、平均を 0 とする正規分布を使うことにすると、以下に示すように面白いことが起こる。

尤度関数 L(data | model) を最大化する最尤推定の枠組みで考えると、

 \begin{eqnarray} L(data | model) &=& P(data | model)\\ &=& \prod_i \frac{1}{\sqrt{2\pi\sigma^2}} e^{-\frac{(y_i - f_i)^2}{2\sigma^2}}\end{eqnarray}

となる。扱いやすくするために対数を取って、対数尤度

 LL(data | model) = \sum_i \left(-\frac{(y_i - f_i)^2}{2\sigma^2}\cdot\log{\frac{1}{\sqrt{2\pi\sigma^2}}}\right)

を最大化するといってもよい。ここで、i によらない定数部分を外に出すと、

 LL(data | model) = -{1\over2\sigma^2}\log{\frac{1}{\sqrt{2\pi\sigma^2}}}\sum_i(y_i - f_i)^2

である。LL は総和記号の中が最大のときに最大になるから、なんと最小二乗法と一致したのである。つまり最小二乗法とは、誤差が正規分布するときの最尤推定に他ならない。

ここで、実験ごとに測定精度が異なる状況を考える。その時は、信頼性の低いデータの重みを軽くして取り扱いたくなる。重み付け最小二乗法では、\sum w_i(f_i - y_i)^2を最小化する。最尤推定の枠組みでは、データごとに誤差のバラつき、すなわち \sigma_i が異なると考えると、{1\over2\sigma_i^2}\log{\frac{1}{\sqrt{2\pi\sigma_i^2}}}が重みに対応している。