最尤推定法 (Maximum Likelyhood や Maximum Likelyhood Estimation と言われ、それぞれ頭文字を取って ML や MLE などとも言われる) は機械学習やコンピュータービジョンなどの分野で良く使われる推定法で、次のような条件付き同時確率を最大化することでパラメータの推定を行います。
$$ \hat{\theta} = \mathop{\mathrm{argmax}}\limits_{\theta} \mathrm{P}(x_1, x_2, \cdots, x_N|\theta) $$
これだけ見て「うん、アレね」と理解できる人はこの記事の対象読者ではなさそうですので、逆にいろいろ教えて下さい。この記事では理論の面から最尤推定法にアタックしてみます。数式成分が多めで、うっとなることもあるかもしれませんが、ゆっくり読んでいきましょう。
※この記事を読むにあたっては、確率論と微分の基礎知識程度は必要です。
尤度
いきなり応用問題から始めると必ず挫折するので、まずは一番簡単な問題設定から始めます。Wikipedia に 最尤推定 のページがありますので、この中で使われている例を参考に話を進めていきましょう。
見た目がまったく同じ 3 枚のコイン A, B, C があります。これらのコインはイカサマで、投げた時に表の出る確率がそれぞれ違います。
- コイン A は 1⁄3 の確率で表が出る
- コイン B は 1⁄2 の確率で表が出る
- コイン C は 2⁄3 の確率で表が出る
この 3 枚のコインを箱の中に入れてシャッフルした後に 1 枚引きます。引いたコインを 10 回投げたら、表が 3 回、裏が 7 回出ました。あなたは A, B, C のどのコインを引いたでしょうか？
この問題設定は極めて単純です。単純すぎて最尤推定法を使わなくても解けるくらい簡単ですが、ここでは敢えて最尤推定法を使って解いてみます。最尤推定法は次の条件付き同時確率を最大化するパラメータ $\hat{\theta}$ を求めることでした。
$$ \hat{\theta} = \mathop{\mathrm{argmax}}\limits_{\theta} \mathrm{P}(x_1, x_2, \cdots, x_N|\ \theta) $$
$x_1$ は 1 回目に投げた時の結果、$x_2$ は 2 回目に投げた時の結果、$x_N$ は N 回目に投げた時の結果になります。コインは全部で 10 回投げることになっているので $N=10$ なのは自明ですよね？これらは、統計学や機械学習の文脈では観測データや学習データとも言われますが、先の問題設定では表が 3 回、裏が 7 回出たということなので、例えばですがこんな風に表せます。(表、裏、みたいな単語だと扱いにくいので、今後は表 = 1、裏 = 0 という風に数値に変換して扱います)
- $x_1 = 0$ (裏)
- $x_2 = 1$ (表)
- $x_3 = 0$ (裏)
- $x_4 = 0$ (裏)
- $x_5 = 0$ (裏)
- $x_6 = 0$ (裏)
- $x_7 = 1$ (表)
- $x_8 = 1$ (表)
- $x_9 = 0$ (裏)
- $x_{10} = 0$ (裏)
ここでちょっとコイントスについて考えます。問題ではコインを 10 回投げていますが、それぞれのコイントスは独立だと考えます。要するに 1 回目に表か裏かどちらが出ようが、2 回目のトスには何の影響も与えないと考えます (3 回目以降も同様) 。まあ普通そうですよね。各試行が独立だと考えると、同時確率は各試行の確率の積に分解できるので、こう書き換えることができます。
$$ \mathrm{P}(x_1, x_2, \cdots, x_{10}|\ \theta) = \prod_i^{10} \mathrm{P}(x_i|\ \theta) $$
じゃあ $\mathrm{P}$ と $\theta$ は何かというと、確率分布とそのパラメータです。確率分布はそれぞれ固有のパラメータを持っていて、そのパラメータが決まれば分布の形が決まります。たとえば…
- 正規分布: 平均 $\mu$ と分散 $\sigma^2$ という 2 つのパラメータを持つ
- t分布: 自由度 $\nu$ というパラメータを持つ
- ベルヌーイ分布: $\lambda$ というパラメータを持つ
式中ではパラメータのことを $\theta$ と書いてはいますが、これを見てわかるようにパラメータは確率分布によって全然違うのでそれを何か 1 つの文字にまとめているというだけで、確率分布によって様々に変わります。
ここでは箱から引いたコインを投げた時にどれくらいの確率で表が出るのか、ということを知りたいので、確率 $\lambda$ で 1 が出て、確率 $1 - \lambda$ で 0 が出る、というベルヌーイ分布を考えます。たとえば $\lambda = 0.4$ だとすると、40% の確率で 1 が、60% の確率で 0 が出ることになります。ベルヌーイ分布はこのように数式で表すことができます。
$$ \mathrm{P}(x\ |\ \lambda) = \lambda^x (1-\lambda)^{1-x} \ \ \ \ (0 \le \lambda \le 1,\ \ x \in \{0,1\}) $$
これを最尤推定法の式に代入すると最終的にはこう書くことができます。
$$
\begin{eqnarray}
\mathrm{P}(x_1, x_2, \cdots, x_{10}|\ \theta) & = & \prod_i^{10} \mathrm{P}(x_i|\ \theta) \
& = & \prod_i^{10} \lambda^{x_i} (1-\lambda)^{1-x_i}
\end{eqnarray}
$$
この式がどういうことを意味しているのかをちょっと考えてみます。この式は 尤度 と呼ばれており、パラメータ $\lambda$ が与えられた時に、その $\lambda$ がどれくらい学習データを尤もらしく説明できているかという指標になります。尤度という概念はとてもわかりづらく、言葉だけでは少し理解しにくいので、実際に値を計算してみましょう。たとえば $\lambda$ を適当に 3 つ決めて、それぞれの $\lambda$ で式の値を求めてみます。
- $\lambda = 0.05$ (5% の確率で表が出るコイン) の場合
$$
\begin{eqnarray}
\prod_i^{10} 0.05^{x_i} (1-0.05)^{1-x_i} & = & 0.95 \cdot 0.05 \cdot 0.95 \cdot 0.95 \cdot 0.95 \cdot 0.95 \cdot 0.05 \cdot 0.05 \cdot 0.95 \cdot 0.95 \
& = & 0.00008729216
\end{eqnarray}
$$
- $\lambda = 0.33$ (33% の確率で表が出るコイン) の場合
$$
\begin{eqnarray}
\prod_i^{10} 0.33^{x_i} (1-0.33)^{1-x_i} & = & 0.67 \cdot 0.33 \cdot 0.67 \cdot 0.67 \cdot 0.67 \cdot 0.67 \cdot 0.33 \cdot 0.33 \cdot 0.67 \cdot 0.67 \
& = & 0.00217803792
\end{eqnarray}
$$
- $\lambda = 0.91$ (91% の確率で表が出るコイン) の場合
$$
\begin{eqnarray}
\prod_i^{10} 0.91^{x_i} (1-0.91)^{1-x_i} & = & 0.09 \cdot 0.91 \cdot 0.09 \cdot 0.09 \cdot 0.09 \cdot 0.09 \cdot 0.91 \cdot 0.91 \cdot 0.09 \cdot 0.09 \
& = & 0.00000003604
\end{eqnarray}
$$
これは、最初に学習データとして得られた $x_1, x_2, \cdots, x_{10}$ が どういうコインから生成されたものか という指標をしめしています。10 回中 3 回だけ表が出るという結果が、$\lambda = 0.91$ のコインから生成されることは稀でしょう。$\lambda = 0.91$ ということは表が出る確率が 91% もあるので、3 回とは言わずもっと表が出るはずです。逆に $\lambda = 0.05$ の場合は、表が出る確率は 5% ということなので、10 回中 3 回も表が出ているという結果はあまり納得のいく結果ではありません。結局は $\lambda = 0.33$ というコインが、一番うまく結果を説明しています。
これは、尤度の値を見てみれば一目瞭然です。3 つの尤度を比べてみると $\lambda = 0.33$ の尤度が最も大きく、続いて $\lambda = 0.05$、そして $\lambda = 0.91$ が最も小さくなっています。要するに尤度の値が大きければ大きいほど、うまくデータを説明してくれる尤もらしいパラメータであると言えます。
この辺りまで理解できるとゴールが見えてきます。最初に示した最尤推定法の式を覚えていますか？これは今考えているコインの問題に対しては、最終的にこのように変形できます。
$$
\begin{eqnarray}
\hat{\theta} & = & \mathop{\mathrm{argmax}}\limits_{\theta} \mathrm{P}(x_1, x_2, \cdots, x_{10}|\ \theta) \
& = & \mathop{\mathrm{argmax}}\limits_{\theta} \prod_i^{10} \mathrm{P}(x_i|\ \theta) \
& = & \mathop{\mathrm{argmax}}\limits_{\lambda} \prod_i^{10} \lambda^{x_i} (1-\lambda)^{1-x_i}
\end{eqnarray}
$$
これは尤度を最大にする $\lambda$ を求める、という意味の式です。尤度が最大になる $\lambda$ がわかれば、観測された学習データ $x_1, x_2, \cdots, x_{10}$ がどういうコインを使って生成されたものなのか、言い換えればどういうコインを投げると $x_1, x_2, \cdots, x_{10}$ という結果が出てくるのか、ということがわかります。
尤度を最大化するパラメータを推定する、これが最尤推定です。
最適化問題
尤度を最大化するという最適化問題なので、あとは最適化問題を解くためのアルゴリズムに従って解いていくだけです。ここでは最大化したい尤度を $L$ として、このように表します。
$$ L(\lambda) = \prod_i^{10} \lambda^{x_i} (1-\lambda)^{1-x_i} $$
ただし、一般的にはこの $L$ をそのまま最大化するのではなく尤度の対数を取った対数尤度 $\log L$ を最大化します。
$$
\begin{eqnarray}
\log L(\lambda) & = & \log \prod_i^{10} \lambda^{x_i} (1-\lambda)^{1-x_i} \
& = & \sum_i^{10} \log \lambda^{x_i} (1-\lambda)^{1-x_i} \
& = & \sum_i^{10} \left( x_i \log \lambda + (1 - x_i) \log (1 - \lambda) \right)
\end{eqnarray}
$$
対数を取る理由は主には、計算を簡単にするためとアンダーフローを防ぐためです。今計算しようとしている尤度は要するに同時確率なので確率の積になっていますが、対数を取ると積を和に変換できます。また、さっき具体的にパラメータを代入して尤度を計算してみてわかったと思いますが、同時確率の計算は値が急激に小さくなっていきます。これも対数を取ることで積が和に変換されるので、値が小さくなることを防げます。
さて、この $\log L(\lambda)$ が最大になるような $\lambda$ を探すことが目的でした。関数の最大値を探す時は微分した結果を = 0 と置いて、微分した変数について解けばよいですよね (関数が単純な凸関数である場合に限られますが) 。
$$
\begin{eqnarray}
\frac{d \log L(\lambda)}{d \lambda} & = & \sum_i^{10} ( \frac{x_i}{\lambda} - \frac{1 - x_i}{1 - \lambda}) \
& = & \frac{1}{\lambda(1 - \lambda)} \sum_i^{10} (x_i - \lambda)
\end{eqnarray}
$$
微分結果が求まったので、これを = 0 と置いて $\lambda$ について解きます。
$$
\begin{eqnarray}
\frac{1}{\lambda(1 - \lambda)} \sum_i^{10} (x_i - \lambda) & = & 0 \
\sum_i^{10} x_i - 10 \lambda & = & 0 \
\lambda & = & \frac{1}{10} \sum_i^{10} x_i
\end{eqnarray}
$$
これで答えがでました。このような式で求めた $\lambda$ が尤度を最大にします。実際に $x_1, x_2, \cdots, x_{10}$ を代入して計算してみると
$$
\begin{eqnarray}
\lambda & = & \frac{1}{10} \sum_i^{10} x_i \
& = & \frac{1}{10} (0 + 1 + 0 + 0 + 0 + 0 + 1 + 1 + 0 + 0) \
& = & \frac{3}{10}
\end{eqnarray}
$$
これは 3⁄10 の確率で表が出るコインであれば尤度が最大になるということで、これに一番近いコインは A になり (A は 1⁄3 の確率で表が出るコイン)、観測された学習データから考えるとあなたが引いたコインは A だろうということが言えます。
気付いた人もいるかもしれませんが、10 回中 3 回表が出るんだったら、そりゃあ 3⁄10 になるのは直感的に当然だろうって思いますよね。これは問題設定が単純ってこともあるかもしれませんが、$x_1, x_2, \cdots, x_{10}$ というデータだけから判断した結果なのでそうなるのも分かる気がしますね。
正規分布のパラメータ推定
最後に正規分布についても例を見て終わりたいと思います。
ベルヌーイ分布のパラメータは $\lambda$ ひとつだけでしたが、正規分布は平均 $\mu$ と分散 $\sigma^2$ の 2 つのパラメータを持っていますので、尤度を最大化するためにそれぞれ $\mu$ と $\sigma^2$ をそれぞれ求めます。今回は確率分布 $\mathrm{P}$ に正規分布を使うので尤度はこうなります。
$$
\begin{eqnarray}
L & = & \mathrm{P}(x_1, x_2, \cdots, x_N|\ \theta) \
& = & \prod_i^N \mathrm{P}(x_i|\ \theta) \
& = & \prod_i^N \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x_i - \mu)^2}{2\sigma^2} \right)
\end{eqnarray}
$$
同じようにこの尤度を $L$ として対数を取ると、このように変形できます。
$$
\begin{eqnarray}
\log L & = & \log \prod_i^N \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x_i - \mu)^2}{2\sigma^2} \right) \
& = & \sum_i^N \log \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x_i - \mu)^2}{2\sigma^2} \right) \
& = & \sum_i^N \left( -\log \sqrt{2 \pi} -\frac{1}{2}\log \sigma^2 - \frac{(x_i - \mu)^2}{2\sigma^2} \right)
\end{eqnarray}
$$
ベルヌーイ分布を考えた時と同様に、この対数尤度関数を $\mu$ と $\sigma^2$ について微分して、それぞれ = 0 と置いて解いていきます。正規分布ではパラメータが 2 つある多変数関数なので偏微分になります。
まずは $\mu$ で偏微分をするところから。
$$
\begin{eqnarray}
\frac{\partial \log L}{\partial \mu} & = & \sum_i^N \left(0 + 0 + \frac{(x_i - \mu)}{\sigma^2} \right) \
& = & \frac{1}{\sigma^2} \sum_i^N (x_i - \mu)
\end{eqnarray}
$$
これを = 0 と置いて $\mu$ について解いていきます。
$$
\begin{eqnarray}
\frac{1}{\sigma^2} \sum_i^N (x_i - \mu) & = & 0 \
\sum_i^N x_i - N\mu & = & 0 \
\mu & = & \frac{1}{N} \sum_i^N x_i
\end{eqnarray}
$$
これでまずは尤度を最大にする $\mu$ が見つかりました。次は $\sigma^2$ について解いていきます。まずは偏微分から。
$$
\begin{eqnarray}
\frac{\partial \log L}{\partial \sigma^2} & = & \sum_i^N \left(0 - \frac{1}{2\sigma^2} + \frac{(x_i - \mu)^2}{2\sigma^4} \right) \
& = & -\frac{N}{2\sigma^2} + \frac{1}{2\sigma^4}\sum_i^{N}(x_i-\mu)^2
\end{eqnarray}
$$
同じようにこれを = 0 と置いて $\sigma^2$ について解きます。
$$
\begin{eqnarray}
-\frac{N}{2\sigma^2} + \frac{1}{2\sigma^4}\sum_i^{N}(x_i-\mu)^2 & = & 0 \
-N\sigma^2 + \sum_i^{N}(x_i-\mu)^2 & = & 0 \
\sigma^2 & = & \frac{1}{N}\sum_i^{N}(x_i-\mu)^2
\end{eqnarray}
$$
これで尤度を最大にする正規分布のパラメータ $\mu$ と $\sigma^2$ がわかりました。
$$
\begin{eqnarray}
\mu & = & \frac{1}{N} \sum_i^N x_i \
\sigma^2 & = & \frac{1}{N}\sum_i^{N}(x_i-\mu)^2
\end{eqnarray}
$$
これは純粋なデータの平均と分散になっています。
まとめ
尤度を最大化してパラメータを推定するための最尤推定法を紹介してきました。学習データからモデルのパラメータ $\theta$ を訓練することが目的ではなく、最終的にやりたいことは新しい未知データがどれだけモデルに適合しているのか推論することです。最尤推定は比較的簡単な手法だと思うので、試しに実装してみると面白いと思います。
