自然科学の統計学4章最尤法を読んだその2

はじめに
記事リンク
4.2 最尤法
続きはまた今度

はじめに

前回の記事はびっくりするくらい誰にも読んでもらえなかった。
僕の目的関数は基本的には承認欲求を最大化するように設計されているので、正直もう続きを書くのはやめてしまおうかとも思ったが、普通に書くのが楽しいのでやっぱり続きを書くことにした。人間はやりたいことをやるべきだし、やりたくないことはやらないべきだ。でも普通にもっと読んでほしい。
あらゆる人間が数学と統計を勉強したい気持ちでいっぱいだということはよく知られた普遍的事実なので、たぶん宣伝が足りなかったのか、冒頭のオタク日記がクソつまんなかったか、冒頭のオタク日記がクソつまんなかったかの3択だと思う。オタク日記を取り下げるつもりはないので、もっとガンガン宣伝していきたい。

記事リンク

第一回/前回
kriver-1.hatenablog.com

次回
kriver-1.hatenablog.com

4.2 最尤法

4.2.1 尤度と対数尤度

前回の記事で、1. データに対して、2. モデルを決めて、3. パラメータを最適化する、という3工程のうち、1と2について説明した。最尤法は3. を担当する手法だ。

データに対してパラメータが「良い」とはどういうことだろうか？即答できる人はもうこの記事を読む必要はない。ブラウザバックして高みを目指す道を邁進してほしい。
ここでは次のような「良さ」を考えてみる。すなわち、そのパラメータからデータがどれだけ出てきやすいか、という指標だ。

例えば、6の面だけがやたらと重いサイコロを振ることを考える。このサイコロは6の面だけがやたらと重いので、反対側である1の目が非常に出やすい。いま、次のようなデータが得られたとしよう。

	1	2	3	4	5	6	全体
出た数	41	10	19	15	12	3	100

データの山 $Y$ は、 $Y=\{(y_1=1),\cdots,(y_{41}=1),(y_{42}=2),\cdots,(y_{51}=2),(y_{52}=3),\cdots,(y_{100}=6)\}$ のように表せるだろう。知りたいのは $y$ の確率分布だ。

普通に考えれば、パラメータの数は6種類、すなわち、 $p_1=(\text{1が出る確率}), p_2=(\text{2が出る確率}), \cdots, p_6=(\text{6が出る確率})$ だ。ただし、合計が1でなければならないので、パラメータの自由度は5次元になってしまっている。
いま、2つの手法を用いて、6つのパラメータを次のように決めてみた。

	1	2	3	4	5	6	全体
出た数	41	10	19	15	12	3	100
手法1	0.50	0.10	0.10	0.10	0.10	0.10	1.00
手法2	0.10	0.10	0.10	0.10	0.10	0.50	1.00

さて、手法1と手法2、どちらがより「それっぽい」だろうか？

このパラメータはどちらも僕が勝手に決めたものだが、手法1は1の目が出やすくなっているのに対して、手法2は6の目が出やすくなっている。データを見ると1の目は6の目の10倍以上も出ているのだから、手法1のほうが明らかに「それっぽい」。このそれっぽさをうまく数式で表してみたい。

そこで、前述した「そのパラメータからデータがどれだけ出てきやすいか」という指標を数式にして、実際に計算してみる。
手法1のパラメータでは、1の目が確率0.5で、2～6の目がそれぞれ確率0.1で出ると主張している。そのパラメータに従うサイコロが仮にあったとして、そのサイコロからデータ $Y$ が得られる確率を求めてみよう。
サイコロを100回振って、各目が $n_1, \cdots, n_6$ だけ出る確率は、多項分布を用いて計算できる。多項分布の実際の式に興味がある人はそんなに多くないとは思うが、せっかくなので途中の計算も載せておくと、
$\displaystyle \begin{eqnarray} Pr(Y; \theta_1) &=& \frac{N!}{n_1! n_2! \cdots n_6!} {p_1}^{n_1} {p_2}^{n_2} \cdots {p_6}^{n_6} \\ &=& \frac{100!}{41! 10! 19! 15! 12! 3!} {0.5}^{41} {0.1}^{10} \cdots {0.1}^{3} \\ &=& \frac{100!}{41! 10! 19! 15! 12! 3!} {0.5}^{41} {0.1}^{59} \\ &\simeq& 7.6 * 10^{-9} \end{eqnarray}$
となる。
同様に、手法2のパラメータに従うサイコロについても計算すると、
$\displaystyle \begin{eqnarray} Pr(Y; \theta_2) &=& \frac{100!}{41! 10! 19! 15! 12! 3!} {0.1}^{41} \cdots {0.1}^{12} {0.5}^{3} \\ &=& \frac{100!}{41! 10! 19! 15! 12! 3!} {0.5}^{3} {0.1}^{97} \\ &\simeq& 2.1 * 10^{-35} \end{eqnarray}$
となる。

途中の煩雑な計算なんかどうだって良いのだが、結局、手法1のパラメータでサイコロを $10^9$ 回振れば1回くらいはデータ $Y$ がそっくりそのまま出てきてくれるが、手法2のパラメータではその ${10}^{26}$ 倍は振らないといけない、ということがわかる。 $10^{26}$ というのは途方もなく大きな数字だ。万億兆京くらいまではわかるかもしれないが、 $10^{26}$ はそのさらに2つ上、100𥝱だ。意味がわからない。

そんなわけなので、この「パラメータからそのデータがどれだけ出てきやすいか」という指標は、パラメータの良さの指標としてまあまあ悪くないものになっているように見える。でかければ良い、小さければ悪い、ということだ。
有用そうなものには名前を付けておきたいので、この指標をこれからは尤度（ゆうど）と呼ぶことにしよう。尤度は英語では likelihood といって（ここでの like は好きではなくそれっぽいという意味だ）、尤もらしさ、とか訳されることもある。君の言うことはもっともだ、とか言うときのもっともである。尤度は英語でlから始まるので $L$ で表すことが多い。

また、尤度は実際の値に興味があるというよりは、でかいか小さいかだけわかればよいという場面が多い。というのも、データが出てくる確率なんてのはデータが増えればそれだけ減っていくわけで、確率そのものに意味なんかないからである。先程の例でも、 $10^{-9}$ や $10^{-35}$ といった、普段なら絶対使わないような桁が出てきていた。
そのかわり、尤度は確率の積で表されることが多い。データが独立に出てくることを考えれば、全体の確率は個々の確率の積になるからである。
積を扱うときは、対数を取ると問題が簡単になることが多い（積が和に分解できるので）。奇跡的にも、対数関数は適用しても大小関係が変わらない、単調増加な関数である。そこで、尤度の対数を取った対数尤度 (log likelihood) が尤度の代わりに用いられることが多い。どうせ尤度自体に意味なんてないので、扱いやすい形に変形してしまうのである。

先程の手法1の尤度を、対数尤度の形で計算してみよう。対数は積を和に分解できるので、
$\displaystyle \begin{eqnarray} \log L(\theta_1; Y) &=& \log Pr(Y; \theta_1) \\ &=& \log \left( \frac{N!}{n_1! n_2! \cdots n_6!} {p_1}^{n_1} {p_2}^{n_2} \cdots {p_6}^{n_6} \right) \\ &=& \sum_{i=1}^{100}\log i - \sum_{j=1}^{6} \sum_{k=1}^{n_j} \log k + \sum_{j=1}^{6} n_j \log p_j \\ &\simeq& -19 \end{eqnarray}$
となる。
計算式がエグいのは変わらないが、積や階乗や累乗が全部ただの足し算と定数倍で表せていることがわかると思う。式中、 $L(\theta; Y)$ は前述の通り尤度 (Likelihood) のLを表している。データが与えられたときの $\theta$ のそれっぽさ、という気持ちを込めて、左に $\theta$ を書いてみた。この辺の順番とかは関数の定義次第なので、記事中で統一されていればどうだって構わない。

手法2についても計算してみると、（数式はまるまる省略するが、） $\log L(\theta_2; Y) \simeq -80$ という値が出てくる。先程も言ったように値そのものに意味はなく、その大小関係が大事なので、まあ手法1のほうがよりそれっぽいんだな、ということがわかる。

4.2.2 最尤推定量と最尤法

さて、我々は対数尤度という強力な指標を得た。また、対数尤度は、パラメータのそれっぽさをなんとなく表す値であるということもわかった。

最尤法の考え方は非常に単純だ。対数尤度が大きいほど良いのだから、対数尤度が最大になるような $\theta$ を考えて、それを最適なパラメータだとしてしまえばいいのである。最も尤度を大きくする、最も尤もらしい $\theta$ を選ぶので、この手法は最尤法と呼ばれている。

$L(\theta)$ を最大にするような $\displaystyle \hat{\theta} = \mathop{\arg\max}_{\theta} L(\theta)$ を、最尤推定量と呼ぶ。推定量 (estimator) というのは、その名の通りパラメータを推定した量のことで、データを入れるとパラメータの予測値が出てくる関数のことでもある。「どんなデータが来ても無視して $\theta=0$ と主張する」というのも、例えば推定量の一種である。最尤推定量は、そんな数多ある推定量の中でも、データに対する対数尤度が最大になるようなパラメータを返す推定量ということだ。

最尤推定量を求めるには、対数尤度をパラメータで微分してあげればいい。関数の最大最小を求めよと言われたら、とりあえず変数で微分するのが常套手段だ。局所解や最小値が出てしまうかもしれないじゃないか、と怒る人もいるかもしれないが、細かいことは気にしなくていい。ほとんどの点は最大でも最小でもないのだから、微分して0になる点がいくつか出てきてくれればそれで御の字だ。どれが最大なのかはその後でゆっくり考えればいい。
なので、対数尤度をパラメータで微分して、それが0になるとした方程式、
$\displaystyle \frac{\partial}{\partial \theta} \log L(\theta) = 0$
のことを、尤度方程式などと呼んだりする。これも大事な式なので、名前をつけたというわけだ。

線形代数を詳しく学んでいない人は、ベクトル量である $\theta$ が微分の分母に来ていることに違和感を覚えるかもしれない。あまり深く考えなくてもよくて、多変数関数 $f(\boldsymbol{x})$ をベクトル $\boldsymbol{x}$ で微分するときは、単に $\boldsymbol{x}$ の各要素 $x_i$ で偏微分して、その結果を元の $\boldsymbol{x}$ と同じように並べてあげればよい。
$\displaystyle \frac{d }{d (x_1, x_2)} f(x_1, x_2) = \left( \frac{\partial }{\partial x_1}f(x_1, x_2), \frac{\partial }{\partial x_2}f(x_1, x_2) \right)$
ということだ。

4.2.3 最尤法の例

せっかくなので、先ほどの違法サイコロのパラメータを最尤推定してみよう。

	1	2	3	4	5	6	全体
出た数	41	10	19	15	12	3	100

パラメータは $p_1, \cdots, p_6$ の6個、制約式として $\sum_i^6 p_i = 1$ があって、対数尤度 $\log L$ は
$\displaystyle \log L(\theta; Y) = \sum_{i=1}^{100}\log i - \sum_{j=1}^{6} \sum_{k=1}^{n_j} \log k + \sum_{j=1}^{6} n_j \log p_j$
で表されるのであった。ここまでは復習だ。

さて、ただの最大最小問題と違って、今回は制約式 $\sum_i^6 p_i = 1$ がある。つまり、そのまま微分しても解は得られない。
こういうときは、2通りの方法がある。1. 制約式を解いてパラメータを減らすか、2. ラグランジュの未定乗数法を使うかである。

4.2.3.1 パラメータを減らす

パラメータを減らすと言っても、特に難しいことはない。 $\sum_i^6 p_i = 1$ を $p_6$ について解いてみるだけである。答えは簡単で、 $p_6 = 1 - p_1 - p_2 - \cdots - p_5$ となる。
$p_1, \cdots, p_5$ を自由に決めれば、 $p_6$ は勝手に決まってしまうということだ。対数尤度は、
$\displaystyle \log L(\theta; Y) = \sum_{i=1}^{100}\log i - \sum_{j=1}^{6} \sum_{k=1}^{n_j} \log k + \sum_{j=1}^{5} n_j \log p_j + n_6 \log (1 - p_1 - \cdots - p_5)$
となる。

これをパラメータ $\theta = (p_1, \cdots, p_5)$ で微分すればいいのだが、先述の通りベクトルで微分するのは各要素で微分して並べるのと同じなので、各要素で微分してそれが0になるという式を考えてしまえばよい。
$\displaystyle \begin{eqnarray} \frac{\partial}{\partial p_j} \log L(\theta; Y) &=& \frac{\partial}{\partial p_j} \sum_{j'=1}^{5} n_{j'} \log p_{j'} + \frac{\partial}{\partial p_j} n_6 \log \left( 1 - \sum_{j'=1}^{5} p_{j'} \right) \\ &=& n_j \frac{\partial}{\partial p_j} \log p_j + n_6 \frac{\partial}{\partial p_j} \log \left( 1-\sum_{j'=1}^5 p_{j'} \right) \\ &=& \frac{n_j}{p_j} + n_6 \frac{-1}{1-\sum_{j'=1}^5 p_{j'}} \\ &=& \frac{n_j}{p_j} - \frac{n_6}{p_6} &=& 0 \\ &\Leftrightarrow&\frac{n_j}{p_j} &=&\frac{n_6}{p_6} \end{eqnarray}$
途中、微分のおかげで $p_j$ を含まない式をばっさりカットできたり、パラメータとして使わないと言っていた $p_6$ を復活させたりもしたが、まあこの式を真面目に追っかけている人はそんなに多くはないと思うので結論だけ述べると、各 $j = 1, 2, \cdots, 5$ について $\frac{n_j}{p_j}=\frac{n_6}{p_6}$ が成り立つ、すなわち各 $j = 1, 2, \cdots, 6$ について度数と確率の比 $\frac{n_j}{p_j}$ が等しくなるということが導かれている。
この比を仮に $k$ とでも置いてあげれば、 $p_j = \frac{n_j}{k}$ となるわけだが、いま $\sum_j p_j = \frac{\sum_j n_j}{k} = \frac{100}{k} = 1$ だったのだから、当然 $k=100$ でなければならない。
すなわち、各パラメータは $p_j = \frac{n_j}{100}$ となったわけである。

結論を並べてみよう。

	1	2	3	4	5	6	全体
出た数	41	10	19	15	12	3	100
最尤推定量	0.41	0.10	0.19	0.15	0.12	0.03	1.00
手法1	0.50	0.10	0.10	0.10	0.10	0.10	1.00
手法2	0.10	0.10	0.10	0.10	0.10	0.50	1.00

うーん、それっぽい。
（最尤推定量が推定量として良いものかどうかの議論をまだしていないので、手放しに良いとは言い切れないのである。）

4.2.3.2 ラグランジュの未定乗数法を使う

蛇足かもしれないが、いまの計算をラグランジュの未定乗数法を使ってやってみよう。
ラグランジュの未定乗数法は非常に便利なワザなので、証明はさておき使えるようになっておくとよい（証明も一般性を欲張り過ぎなければ非常に単純だ）。こういう死ぬほど便利な☆5スキルを、ガチャを引かなくても全人類が習得できるのが現実世界のいいところだ。

ラグランジュの未定乗数法では、まず制約式を $f(x) = 0$ の形に直す。今回の制約式は $\sum_j p_j = 1$ なので、 $f(\theta) = \sum_j p_j - 1 = 0$ となる。
つぎに、いまの制約式に未定乗数 $\lambda$ を掛けて、目的関数（最大・最小を求めたい関数）に足す。今回の目的関数は対数尤度 $\log L(\theta; Y)$ なので、変形後の目的関数は
$\displaystyle F(\theta, \lambda) = \log L(\theta; Y) + \lambda \left( \sum_j^6 p_j - 1 \right)$
のようになるだろう。
あとはこれを、制約式のない普通の最大・最小問題だと思って微分すればよい。これでうまくいくよ、というのがラグランジュの未定乗数法である。

魔法のような話かもしれないが、信じられないなら実際にやってみればよい。先ほどと同様に $F$ を各パラメータで微分する（ここで、新しい目的関数 $F$ にはパラメータ $\lambda$ が追加されていることを忘れてはならない）。

$\displaystyle \begin{eqnarray} \frac{\partial}{\partial p_j} \log F(\theta, \lambda) &=& \frac{\partial}{\partial p_j} \sum_{j'=1}^{6} n_{j'} \log p_{j'} +\frac{\partial}{\partial p_j} \lambda \left( \sum_{j'}^6 p_{j'} - 1 \right) \\ &=& n_j \frac{\partial}{\partial p_j} \log p_j + \lambda \\ &=& \frac{n_j}{p_j} + \lambda &=& 0 \\ &\Leftrightarrow& \frac{n_j}{p_j} &=& -\lambda \end{eqnarray}$

$\displaystyle \begin{eqnarray} \frac{\partial}{\partial \lambda} \log F(\theta, \lambda) &=& \frac{\partial}{\partial \lambda} \lambda \left( \sum_{j}^6 p_{j} - 1 \right) \\ &=& \left( \sum_{j}^6 p_{j} - 1 \right) &=& 0 \\ &\Leftrightarrow& \sum_{j}^6 p_{j} &=& 1 \end{eqnarray}$

これは先ほどの、パラメータを減らすやり方において、 $k = -\lambda$ とした場合と全く同じである。全く同じであるから、当然結果も全く同じであって、各パラメータは $p_j = \frac{n_j}{100}$ と求まる。

どちらのやり方でやったって構わないが、パラメータを減らすのが常に簡単に実行できるとは限らないのに対して、ラグランジュの未定乗数法は関数を追加して微分するだけなので基本的にいつでも使える。☆5スキルとはそういうものなので、ラグランジュ先輩に感謝の念を示しつつ、おとなしくぶっ壊れスキルの強さを享受するのがいいんじゃないだろうか。

続きはまた今度

今回は最尤法というものがどういうものなのかについて説明し、その具体的な計算例を述べた。次回以降は最尤法の良さの説明に必要な「推定量の良さ」の話をいくつかできたらいいなと思っている。
正直この章は次項と次次項がメインで、ここまでは前座でしかないので、前座に対して既に2記事も使っていることに不安がないわけではないが、失踪しない程度にがんばっていきたい。

僕が失踪するかどうかはみんなの反応にかかっている。下にある黄色い星を押すか、その右下にある青い鳥を押してほしい。それが僕の幸せなんだ。よろしくお願いします。

第一回/前回
kriver-1.hatenablog.com

次回
kriver-1.hatenablog.com