最尤法解説の第4回目です。今までの記事（第1回、第2回、第3回）では、最尤法のアイデアとポアソン分布についての解説が終わりました。
軽くまとめておきますと、
という感じでした。
この記事では、今まで確認してきたアイデアを基に、データから確率密度関数の具体的なパラメータを求める方法を考えてみます。そして、その方法が最尤法と呼ばれるものになります。
ところで、お忘れの方もいらっしゃるかもしれませんし、理解していないと最尤法が結局何をやっているのか分からなくなってしまいかねないので、もう一度最尤法のアイデアをまとめておきます。
そもそも達成したいことは、手元にあるデータを利用して母集団の確率密度関数を求めることでした。当然ながら、データから求めた確率密度関数と母集団の確率密度関数が一致してくれれば「目標を達成できた」と言えそうです。
ですが言うは易く行うは難しです。というのも、データから確率密度関数を求めた後、それが母集団の確率密度関数と一致していると言えるためには、そもそも母集団の確率密度関数が分かっている必要があります。ですが、実際には母集団の確率密度関数を知ることは基本的に不可能です。これが”行うは難し”の理由ですね。
「目標を達成できた」と言えるための基準を母集団の確率密度関数と設定している限り、どうあがいても目標を達成できないということですな。
そこで別の考え方が必要になります。その”別の考え方”の一つが、「手元のデータを生成するような確率密度関数こそが母集団の確率密度関数だ」という考え方になります。
その理屈も軽く確認しておきましょう。まず次の2つの考え方がまず成り立っていると考えられます。「現実に実現するデータは、最も高い確率を与えられたデータである」「今手元にデータが存在できてしまっている」という2つの考え方です。
これら2つの考え方が正しければ、「あるデータが今手元に存在できてしまっている（生成されてしまっている）以上、母集団の確率密度関数では、今手元にあるデータを生成する確率が最も高くなっていたのだろう」と考えるのが妥当だと言えそうです。ここから、「手元にあるデータを生成するような確率密度関数こそが母集団の確率密度関数だ」が導けることになります。
このような論理構造で、「手元のデータを生成する確率」という新たな基準を（母集団の確率密度関数ではない基準を）設けることが出来ました。
そこで、今回はこのアイデアに基づいて最尤法の具体的な方法論を確認していきます。
※毎度のことながら注意を。この記事シリーズでは、本来なら「確率質量関数」と表記すべきところも「確率密度関数」と表記しています（例えば、ポアソン分布は確率質量関数と呼ぶべきですが、この記事シリーズでは一貫して確率密度関数と呼んでいます）。それは、最尤法は確率”質量”関数でも確率”密度”関数でも、統一的に扱うことが出来るからです。読んでいて違和感を感じる方もいらっしゃるかと思いますが、ご了承ください。
データを生成する確率はどう考えればいい？
今まで「手元のデータを生成する確率」とさんざん言ってきました。しかし、一言で”データ”といっても、その実態は様々でしょう。単一の数値（例えば、誰かの身長や誰かの最高血圧など）を表すこともあれば、単一の数値の集合体（データセット）を表すこともあるでしょう。場合によっては、単一の数値やその集合体の総称として”データ”という単語を使うこともあるかもしれません。
“データ”に関する補足（この補足は読んでおくことをおすすめします）
記録されてさえいれば、それはどのような性質のものであれ”データ”と呼べてしまう以上、感情の変化や人間の動作などもデータと呼べます。その意味では、ここで「単一の数値」「単一の数値の集合体」といったように、”データ”があたかも数値のみに限定されるかのように書くのは本来不適切です。
しかし、最尤法はその性質上、数値として記録されたデータ（あるいは、何かしらの方法で数値化したデータ）にしか適用できません。そのため、個人的には”データ”が数値以外の様々なものを表し得るとして煩雑な一般論を展開する必要はないと思います。
この記事シリーズの目的は最尤法を解説することです。そのため、”データ”が数値に限定されないという前提の下に煩雑な一般論を展開して話の本筋を分かりにくくしてしまうより、”データ”が数値に限定されるという前提をおいて話の本筋を分かりやすくした方が、より理に適っていると判断しました。
そのような考えから、この記事シリーズでは以降、”データ”と言えば数値で記録されたデータに限定されるという前提を設けさせていただきます。
しかし、現実には数値データだけでなく、動画や画像、文章などといったように、（数値として表現できたとしても）本質的には数値でないようなデータも存在することはご注意ください。
異論もあるかとは思いますが、この記事シリーズではそのような立場で議論を進めていきます。ある程度の厳密性を犠牲に分かりやすさを重視しようとしたために、このような表現になっているのだとご理解ください。
“データ”についての補足 ここまで
一般的にデータという単語が表すものとしては「単一の数値」「単一の数値の寄せ集め」「数値そのものと、その寄せ集めの総称」という3つの意味で使われるでしょう。今までこの記事シリーズでは”データ”という単語を無批判に使ってきましたが、今までは暗に3つ目の意味「数値やその集合の総称」で”データ”という単語を使ってきたということになります。
しかし、「データの生成される確率」を考える上では、その”データ”という単語が「単一の数値、もしくは数値の集合を表す」というように曖昧さを孕んでしまっていては詳細な議論が出来ません。なので、そこで使われている”データ”という単語が「単一の数値のみ」を表すのか「数値の集合のみ」を表すのか、どちらを表しているのかを考えます（考え抜いた結果、やはりどちらも表すとするしかないという結論も一つの重要な結論となります）。
そもそも統計調査である以上、いくつかの数値が必要です。そしてその”いくつかの数値”は、最高血圧と最低血圧というように、単一のデータが複数種類あるという意味ではなく、最高血圧であれば最高血圧という同一種類の単一のデータがいくつもあるという意味であるはずです。
そう考えれば、今まで無批判的に使ってきた「今手元にあるデータが得られる確率」は、「今手元にあるような数値の集合が得られる確率」と考えるべきです。つまり、同一種類の数値が大量にあるという状況下で、そのような数値の集合を生成する確率をどうにか評価する方法を考えなければいけないということになります。
ということで、「今手元にあるデータセット（データではなく）が得られる確率」をどのようにすれば評価できるかを考えていきます。以降、データは「単一の数値」を、データセットは「単一の数値の集合」をそれぞれ表すこととします。
「データ（単一の数値）」ではなく「データセット（数値の集合）」を考えるべきだという話になったばかりですが、ここでまた「データ」の得られる確率を考えます。今は確率密度関数の数式は分かっているという前提があります。
なので、後はパラメータさえ分かれば、その「データ」の得られる確率を計算できるということになります。ここではポアソン分布を考えていますから、パラメータとは平均値\( \lambda \)のことになります。
ここでとりあえず、パラメータを\( \lambda \)、得られたデータ（数値）を\( x \)として、その\( x \)が得られる確率をポアソン分布\( p(x, \lambda) \)と表すことにします。この\( x \)はポアソン分布に従うとしているので、次の(1)式のように表すことが出来ます。
$$ p( x, \lambda ) = \frac{ \lambda ^ x exp ( – \lambda ) } { x! } ・・・ (1) $$
現状ではまだパラメータが分からないので、この確率を具体的に計算することはできませんが、\( x=1 \)となる確率は次の(1′)式のように表せます。
$$ p( 1, \lambda ) = \frac{ \lambda ^ 1 exp ( – \lambda ) } { 1! } = \lambda exp(-\lambda) ・・・ (1′) $$
次に「データ」ではなく「データセット」を考えてみます。データセットを\( X \)（大文字のX）で表し、データセットは\( N \)個のデータで構成されていて、各データは\( x_i \)と表すとします。すると、そのデータセットを得られる確率は次のように表されるはずです。
$$
\begin{equation}
\begin{split}
p(X, \lambda) &= p(x_1) \times p(x_2) \times \ldots \times p(x_i) \times \ldots \times p(x_N) \\
&= \prod_{i=1}^{N} p(x_i) \\
&= \prod_{i=1}^{N} \frac{ \lambda ^ {x_i} exp ( – \lambda ) } { x_i! } ・・・(2) \\
\end{split}
\end{equation}
$$
ぱっと見は複雑そうに見えるのでビビるかもしれませんが、各データの得られる確率を掛け算しただけです。要するに\( x_1 \)が得られて、かつ\( x_2 \)が得られて、・・・かつ\( x_N \)が得られたという同時確率を計算してるわけですな。
途中の\( \prod_{i = 1}^{N} \)という記号は見慣れない人もいらっしゃるでしょうが、これは単に添字iが1番のものから順にN番までのものすべてを掛け合わせるということを表しています。総和記号の積版だと考えていただければ分かりやすいかと。
先ほど「データ」に対する確率は、まだ具体的には計算できないと説明しましたが、やはり同様に、「データセット」に対する確率も計算できません。ですが、先ほどの「データ」に対する確率の式（1）を見てください。その中で分かっていない数値はパラメータ\( \lambda \)と具体的な確率\( p(x, \lambda) \)です。
そこで、試しにパラメータ\( \lambda \)に具体的な数値を色々と代入してみて、確率\( p(x, \lambda) \)がどのように変化するかを見てみましょう。
馬に蹴られて死んだ兵士の数
この記事シリーズでは、具体的な数値としてはボルトキーヴィチによる馬に蹴られて死んだ兵士の数を使うとしていました。まずはそこで報告されている数値を引用してみます。
0人：109件
1人：65件
2人：22件
3人：3件
4人：1件
5人以上：0件
ここで、前回記事の最後の方の注意と同じような注意が必要です。これはあくまでもデータの集計結果であって、データでもデータセットでもないということです（この記事のデータとデータセットの定義によれば）。
馬に蹴られて死んだ兵士の数を調べた場合、データは「馬に蹴られて慎が兵士の数は、今年は○○回」というものになって、そのデータを集めた「今年は○○回、今年は△△回、・・・」というものがデータセットとなります。
そして、そのデータセットから同じ回数のデータをカウントした（集計した）結果が上の数値というわけです。そして、このカウントした結果がポアソン分布に従います。
先ほどの確率計算の式(1)と(2)はどちらも「データ」や「データセット」を対象にしていました。そのため、「集計結果」であるこの数値は直接(1)や(2)に代入することはできません。そのため、集計元となったデータを復元する必要があります。\( x_1 = 109 \)、\( x_2 = 65 \)、・・・となるわけではないという点にご注意ください。
具体的に復元してみましょう。
まず、0人が109件ということは、便宜的に\( x_1 \)から\( x_{109} \)までがすべて0と考えられます。なぜ”便宜的に”なのかと言えば、元のデータセットでは本当に1番目のデータから109番目までのデータがすべて0というようになっていたとは限らないからです。
もしかしたら、集計前のデータでは\( x_1 \)（最初のデータ）は0であっても、その次の\( x_2 \)は1など0以外の数だったかもしれません。つまり、上の集計結果からは、順番までは復元できないということですね（だから、ここに復元という言葉を使うのは本当は不適切かもしれません）。なので、便宜的に（順番は違うかもしれないけど）\( x_1 \)から\( x_{109} \)までは0だったということにしています。
そして同様に\( x_{110} \)から\( x_{174} \)までは1、・・・としていきます。すると、結局元のデータセット（ただし、順番という時系列的な観点からすれば正しくない）としては200個のデータが復元できます。この復元できたデータ200個を先ほどの(2)に代入すれば、データセットの得られる確率が分かるはずです。
具体的に(1)に代入してみましょう。すると、次の(3)式のようになります。
$$
\begin{equation}
\begin{split}
p( X, \lambda ) &= \prod_{i=1}^{N} p(x_i, \lambda) \\
&= \prod_{i=1}^{N} \frac{ \lambda ^ {x_i} exp ( – \lambda ) } { x_i! } \\
&= \prod_{i=1}^{109} \frac{ \lambda ^ {0} exp ( – \lambda ) } { 0! } \\
&\hspace{10pt} \times \prod_{i=110}^{174} \frac{ \lambda ^ {1} exp ( – \lambda ) } { 1! } \\
&\hspace{10pt} \times \prod_{i=175}^{196} \frac{ \lambda ^ {2} exp ( – \lambda ) } { 2! } \\
&\hspace{10pt} \times \prod_{i=197}^{199} \frac{ \lambda ^ {3} exp ( – \lambda ) } { 3! } \\
&\hspace{10pt} \times \prod_{i=200}^{200} \frac{ \lambda ^ {4} exp ( – \lambda ) } { 4! } ・・・(3)
\end{split}
\end{equation}
$$
となります。何をしているかと言えば、先ほど復元したデータを順番通りに代入しています。1番から109番までは\( x_i = 0 \)だったことは確認しました。なので、1番から109番までのデータに対する確率はすべて\( x_i = 0 \)として計算できることになります。110番以降についても同様です。
この(3)式を使って、具体的に色々なパラメータを代入して計算してみます。
確率計算のシミュレーション
まず仮に、パラメータが\( \lambda = 0.1 \)だったとします。その場合、(3)式は次のように計算できます。長いですが、ただ先ほどの(3)式で\( \lambda \)となっていた部分を0.1に変えただけです。丁寧に書き下したせいで、見た目には複雑そうな式になっているだけです。
$$
\begin{equation}
\begin{split}
p( X, 0.1 ) &= \prod_{i=1}^{N} p(x_i, 0.1) \\
&= \prod_{i=1}^{N} \frac{ 0.1 ^ {x_i} exp ( – 0.1 ) } { x_i! } \\
&= \prod_{i=1}^{109} \frac{ 0.1 ^ {0} exp ( – 0.1 ) } { 0! } \\
&\hspace{10pt} \times \prod_{i=110}^{174} \frac{ 0.1 ^ {1} exp ( – 0.1 ) } { 1! } \\
&\hspace{10pt} \times \prod_{i=174}^{196} \frac{ 0.1 ^ {2} exp ( – 0.1 ) } { 2! } \\
&\hspace{10pt} \times \prod_{i=196}^{199} \frac{ 0.1 ^ {3} exp ( – 0.1 ) } { 3! } \\
&\hspace{10pt} \times \prod_{i=200}^{200} \frac{ 0.1 ^ {4} exp ( – 0.1 ) } { 4! } \\
&= 9.48 \times 10^{-142}・・・(3)
\end{split}
\end{equation}
$$
計算機でも使って計算してやると、\( \lambda = 0.1 \)だった場合にボルトキーヴィチが報告したようなデータになる確率は\( p(X, \lambda ) = 9.48 \times 10^{-142} \)ということが分かります。ちなみに、こちらのExcelファイルで計算できますんで、よろしければ参考にしてくださいませ～（一応グラフも付けておきました）。これはかなり小さい確率に見えますが、200個の確率を掛け合わせているので当然の結果です。
長ったらしくなるので、もう計算式は書きませんが、\( \lambda = 0.5 \)だった場合は\( p( X, \lambda = 3.22 \times 10^{-90} ) \)となり、\( \lambda = 0.1 \)だったときよりも大きくなっています。
さらに\( \lambda \)を大きくして、\( \lambda = 1.0 \)だった場合を計算してみると\( p( X, \lambda ) = 6.36 \times 10^{-98} \)となり、\( \lambda = 0.5 \)だったときよりも小さくなっています。
このように、パラメータ\( \lambda \)を小さい側から順に大きくしていくと、ある所までは大きくなり、それ以降は小さくなっていくという傾向が見て取れます。
\( p( X, \lambda ) \)とは、「今手元にあるデータセットが得られる確率」だったわけですから、この数値が一番大きくなるような\( \lambda \)を選べば、この記事シリーズを通しての悲願であった「母集団の確率密度関数を求める」が達成できそうです。
そして、そのような方法こそが最尤法になります。
3つの場合でのシミュレーションの結果から、偶然にも最大となる\( p( X, \lambda ) \)に目星をつけることが出来そうです。この\( p(X, \lambda) \)は\( 0.5 \lt \lambda \lt 1.0 \)辺りで最大値となりそうです。
しかし、このようなあてずっぽうな方法ではとても非効率的だということは簡単に分かると思います。
例えば、いつでも\( \lambda = 0.1 \)から順に計算していけばいいとは限りません。場合によっては1や10といったもっと大きい数値から始める必要があるかもしれません。そのように考えれば、\( p(X, \lambda) \)をもっと効率的に求める方法がほしいと思うものでしょう。
さらに、今回はポアソン分布に限定していますが、他の分布についても統一的に扱えるようにしたいです。
今の計算をさらに洗練した方法が最尤法になります。ただ、その前にはまだ尤度関数という概念を説明する必要があります。
まとめ
しまった。詳しく説明していたら、肝心の尤度関数にたどり着かなかった。とは言え、長くなってきたので、ここらへんでいったんまとめを。結局のところ、次のことが理解できていれば十分です。
- 今まで言ってきた「今手元にあるデータを得られる確率」は、本来「今手元にあるデータセットを得られる確率」と解釈すべきだ
- 今手元にあるデータセットを得られる確率\( p(X, \lambda) \)は、そのデータセット中に含まれるデータ数を\( N \)、確率密度関数を\( p(x_i), \lambda \)として\( \prod_{ i = 1 }^{ N } p( x_i, \lambda ) \)と表される
- 集計結果（今回のボルトキーヴィチのような）を利用するときには、元のデータをある程度復元する必要がある
次回こそは尤度関数について説明できると思います。
P.S. どうでもいいことですが、アイデアを形にするのって楽しいですよね。