はじめに
今回は具体的な手法の詳しい解説やコードなどは一切出てきません。 ある程度の数式（確率分布や微分積分、線形代数）に抵抗の無い方が、最尤推定とベイズ推定はどういうものであるのかを初めて学ぶのに無理の無い内容になっていると思われます。
データと確率分布
まずは基本事項としてデータと確率分布について説明します。
確率論でデータを扱う場合には、データはとある真の確率分布 $\phi(\cdot)$ から生じているものであると考えます。具体的には、今着目しているデータ $x$ は確率変数 $X$ が $\phi(X)$ に従っており、その実現値として $X = x$ と定まったことによって得られていると考えます。このことを
$$ x \sim \phi(X) $$
と表現します。更に確率変数が複数あるケースを考えましょう。ここでは $X, Y$ の2つを考慮します。複数の確率変数がある場合はそれらをまとめた場合の確率分布、"同時分布"を考えることになります。具体的には
$$ X, Y \sim \phi(X, Y) $$
という表現をするということです。ここで $X, Y$ には真の実態として関係があっても無くても良いです。もしも $X$ と $Y$ が互いに全く無関係である場合には同時分布として考えたものが $\phi(X, Y) = \phi _ X(X)\phi _ Y (Y)$ と分解されます。逆にこのように分解できるときには確率変数 $X$ と $Y$ が独立であると表現します。
さてここで、確率変数列 $X _ 1, \cdots, X _ N $ を考えてみます。この確率変数 $X _ i$ の添字 $i$ はデータを取得する順番でつけられているとしましょう。
つまり、今からデータを $N$ 個だけ計測し、$X _ i = x _ i$ として計測結果を埋めていく作業をするということになります。すると、
$$ X _ 1, \cdots, X _ N \sim \phi(X _ 1, \cdots, X _ N) $$
のような真の同時分布から生成されるものとして扱うことになります。今 $X _ i$ と $X _ j$ が互いに独立であるとすれば、既に見たようにそれぞれの確率変数の分布の積で分解できるという話になるわけです。
具体的な例を考えます。ある同じ性質であると考えられる確率変数を複数回観測する場合、例えば、サイコロを $N$ 回振ったときの出目を記録していくことを考えた場合、$i$ 回目の出目の確率変数を $X _ i$ とすると、各回独立であろうというのは何となく想像がつくところでしょう。
しかしこれは単に人間側がそう思っているという話だけであり、真の実態を我々は知りません。 実はデータの生成規則（ここではサイコロの出目の規則）は時刻に応じて性質を変えてしまっているのではないか？ということは念頭に置いておくことは損にならないはずです。例えば、サイコロの重量が凄まじく、投げる度にサイコロの角が消し飛びが削れているということがあるとすれば、投げる度にその確率分布は変化していることが考えられます。すると、必ずしも $X _ 1 \cdots, X _ N$ の同時分布を個々に分解はできず、あくまで同時分布 $\phi(X _ 1, \cdots, X _ N)$ という形式でしか表現のしようがないかもしれません。
いずれにせよ、データを採取したときにはそれを確率変数の実現値としてみなし、背後には真の分布$\phi(\cdot)$ があると考えます。
データと確率モデル
さて、真の分布から生成されたデータ $D = \{x _ 1, \cdots, x _ N\}$ について考えます。真の分布について知る由もない我々は、データを可視化してみたり、あるいは物理的な性質を加味したりして、代わりに何らかの自分たちが取り回しやすい確率分布 $p(\cdot)$ をこしらえることにします。このようにこしらえた確率分布 $p(\cdot)$ のことを確率モデルと呼びます。
特に、確率分布として性質が調べられている固有名詞がつくような分布、例えば正規分布やポアソン分布などを用いると、後々の推定がとてもシステマチックに実施できます。推定とは、データの背後にある真の分布 $\phi(\cdot)$ になるべく近くなるように確率モデル $p(\cdot)$ を調整する作業です。もしも正規分布を選んだ場合、その形状は平均パラメータ $\mu$ と分散パラメータ $\sigma ^ 2$ の2つだけで決定されることになり、無数にある確率分布の形状の中からまともに使える形状を決める作業が、たった2つのパラメータを調整する作業に簡略化できるわけです。
ここで、当然真の分布が正規分布で表現できるとは限らないことに注意しましょう。それでも、正規分布で表現できる範囲の中で分析や予測ができるのであれば、当てずっぽうにデータを解釈したり予測したりするよりは便利な道具として確率モデルを使えるはずです。そういうわけで、確率モデルはデータを自分たちの応用目的に応じて適宜設定されるものであることは念頭においておくべきです。
推定方法
確率モデルの形状を調整する具体的な方法の例を紹介します。
まず、前提として、あるパラメータ $\theta$ を持つ確率モデル $p(\cdot \mid \theta)$ を勝手に設定し、真の分布から生成された $D$ を参考に確率モデルを調整していくことを考えています。また、今回は都合良く、データ1つ1つが独立同分布から生成されるような条件の中で確率モデルを調整していきます。
確率変数 $X _ i$ と $X _ j$ が独立同分布に従うというのは、データの生成規則である分布 $p(\cdot)$ を共有しているということです。すなわち、生成の様子を
$$ \begin{align} X _ i &\sim p(X) \\ X _ j &\sim p(X) \end{align} $$
と考えているということです（同じ $p(X)$ から繰り返し生成されているということ）。
すなわち確率変数の実現値であるデータ $D = \{x _ 1, \cdots, x _ N \}$ に関して、実現値として $X _ i = x _ i$ が得られる確率が
$$ p(X _ 1 = x _ 1, \cdots, X _ N = x _ N) = \prod _ {i = 1} ^ N p(X _ i = x _ i) $$
と計算できることになります。
最尤推定法
確率モデル $p(\cdot \mid \theta)$ をこしらえて、その確率モデルから乱数を発生させてみることを考えます。すなわち疑似データ $a _ i$ を大量に発生させてみるのです。
$$ a _ i \sim p(\cdot \mid \theta) $$
ということを実験してやってみたときに、都合よく手元のデータの実現値 $D = \{x _ 1, \cdots, x _ N \}$ と一致する確率はどれくらいでしょうか。これを実際にデータと一致するような結果が出るまで繰り返して確率を求めても良いかもしれませんが、それは途方もない作業でしょう（そもそも勝手に作った確率モデルに従う乱数生成器なんて作れるとは限らん）。
数式を使って確率モデル $p(\cdot \mid \theta)$ が $D$ を生成する確率を記述して見る方針を取ります。独立同分布を前提にすると
$$ p(D \mid \theta) = \prod _ {i=1} ^ {N} p(x _ i \mid \theta) $$
と表せます。すなわち、具体的に $D$ を生成する確率は単に確率モデル $p(\cdot \mid \theta)$ に $x _ i$ を代入して計算してみれば良いというだけなのです。そうすると、$\theta$ の値によって、この確率モデルがデータを実際に生成するような確率が変わることが想像できます。
今、我々は真の分布 $\phi (\cdot)$ から生成されたデータを目の前に、確率モデル $p(\cdot \mid \theta)$ を勝手に作ってこれを調整しようとしているのでした。そこで、"調整され終えた確率モデルが、既に得られているデータを生成する確率は高くて然るべきである"という考えに基づいているのが最尤推定です。
確率モデル $p(\cdot \mid \theta)$ が既に得られているデータを生成する確率を「尤度」と呼び、尤度が最も高くなるように $\theta$ を調整するのが最尤推定であると言えます。尤度は
$$ L _ D (\theta) = p(D \mid \theta) = \prod _ {i=1} ^ {N} p(x _ i \mid \theta) $$
と表記します。従って最尤推定とは
$$ \hat \theta = {\rm argmax} _ \theta L _ D (\theta) $$
を求める手続きにより実行されます。この最適化問題は確率モデルが単純であれば解析的に求まります。一方で、多くの応用モデルでは数式で求めることが困難であるため、応用上は当てずっぽうの $\theta ^ {(0)}$ から初めて以下の更新規則によってパラメータを調整していきます。
$$ \theta ^ {(i)} \leftarrow \theta ^ {(i-1)} - \alpha \frac{\partial L _ D (\theta)}{ \partial \theta} $$
ここで $\alpha$ を学習率と呼び、上記の更新則を勾配法と呼びます。更新が十分に行われなくなったところで更新を打ち切り、そのときの値を $\hat \theta$ として利用するという形になります（これが本当に最尤推定解である保証は、尤度関数が凸の場合にしか得られないことには注意が必要）。
こうして真の分布は分からないが、
$$ p(X \mid \hat \theta ) $$
という確率分布を作れば、まがい物なりには良い分布が作れたというわけです。
ベイズ法
尤度関数
$$ L _ D (\theta) = p(D \mid \theta) = \prod _ {i=1} ^ {N} p(x _ i \mid \theta) $$
に関して、一度、データが元々確率変数であったことを思い出すことにしましょう。再度データ $D$ を取得し直すことを考えたら、同じ推定方法を用いても推定される $\theta$ は変わることが想定できます。すなわち、$\theta$ はデータ $D$ が確率変数から実現値に変わる過程まで遡ると確率変数であると言えます。そこでいま一度 $\theta$ を確率変数として扱い、同様に $N$ 個の確率変数の実現値であるデータ $D$ も確率変数 $D'$ として取り回すと、
$$ p(D', \theta) $$
という同時分布を考えることできます。同時分布は
$$p(X, Y) = p(X \mid Y)p(Y) = p(Y\mid X)p(X)$$
と
$$\int _ Y p(X , Y) \mathrm d Y = p(X)$$
という関係が常に成り立ちます。これを考慮すると
$$ p(D') = \int _ \theta p(D', \theta) \mathrm d \theta = \int _ \theta p(D' \mid \theta)p(\theta) \mathrm d \theta $$
という、確率モデル $p(\cdot \mid \theta)$ から $\theta$ を積分消去したデータの確率分布を考えることができます。ここで確率変数 $D'$ に対して実現値 $D$ を代入すれば、尤度関数の如く
$$ p(D) = \int _ \theta p(D, \theta) \mathrm d \theta = \int _ \theta \prod _ {i=1} ^ {N} p(x _ i \mid \theta)p(\theta) \mathrm d \theta $$
とデータを生成する確率を計算することができます。この $p(D)$ のことを周辺尤度と呼びます（実は積分消去をすることを周辺化と呼ぶ。なので、同時分布からパラメータ $\theta $ を積分消去したあとのデータの生成確率を表す関数を周辺尤度と呼ぶのである）。
今、将来得られる1つの確率変数 $X$ に着目すると、その $X$ が得られる確率は、上記までのパラメータの周辺化の流れをなぞると
$$ p(X) = \int _ \theta p(X, \theta) \mathrm d \theta = \int _ \theta p(X \mid \theta)p(\theta) \mathrm d \theta $$
と表せそうです。これを"ベイズ（事前）予測分布" と言います。この予測分布を上手に決める事ができれば嬉しいというのがベイズ法のモチベーションです。この予測分布には既に $\theta$ という調整役のパラメータが表面上消えてしまっています。その代わりに、最右辺に、$p(\theta)$ というパラメータの分布がいる状態です。この分布のことを事前分布と呼びます。
ベイズ法では、確率モデル $p(\cdot \mid \theta)$ と事前分布 $p(\theta)$ の2つの組を合わせて（ベイズ）確率モデルと呼びます。最尤推定では $\theta$ を調整していました。一方でベイズ推定では、当てずっぽうの事前分布 $p(\theta)$ をデータ $D$ に基づいて調整します。$\theta$ は確率変数なのでしたから、再度確率変数である $D'$ を持ち出せば
$$ p(D', \theta) = p(\theta \mid D')p(D') = p(D' \mid \theta)p(\theta) $$
という同時分布の関係式が利用でき、
$$ p(\theta \mid D') = \frac{p(D' \mid \theta)p(\theta)}{p(D')} $$
という式変形が実施できます。ここで $D'$ に具体的な実現値である $D$ を代入してやると
$$ \begin{align} p(\theta \mid D) &= \frac{p(D \mid \theta)p(\theta)}{p(D)} \\ & = \frac{\prod _ {i=1} ^ {N} p(x _ i \mid \theta)p(\theta)}{\int _ \theta \prod _ {i=1} ^ {N} p(x _ i \mid \theta)p(\theta) \mathrm d \theta} \end{align} $$
という式に落ち着き、確率変数 $D'$ が $D$ として実現した上での $\theta$ の確率分布である "条件付き分布" $p(\theta \mid D)$ が得られたことになります。この式をちょっとだけ見返しましょう。下記のようにデータに依存している因子とそうでない因子を分けて書いてみます。
$$ p(\theta \mid D') = p(\theta) \times \frac{\prod _ {i=1} ^ {N} p(x _ i \mid \theta)}{\int _ \theta \prod _ {i=1} ^ {N} p(x _ i \mid \theta)p(\theta) \mathrm d \theta} $$
と、事前分布 $p(\theta)$ に対してデータに依存する項 "尤度 / 周辺尤度" を掛けて修正しているように見えます。（これは単に同時分布の関係式そのものなのだが、この様子を特別視してか）上記の式をベイズの定理と呼びます。またこの更新規則をベイズ更新と呼びます。
また、データに条件付けられた $\theta$ の分布が事前分布をデータによって更新して得られていることからか、$p(\theta \mid D)$ を事後分布と呼びます。ベイズ（事前）予測分布
$$ p(X) = \int _ \theta p(X \mid \theta)p(\theta) \mathrm d \theta $$
の $p(\theta)$ を事後分布 $p(\theta \mid D)$ で置き換えたものを
$$ p(X \mid D) = \int _ \theta p(X \mid \theta)p(\theta \mid D) \mathrm d \theta $$
と書き、ベイズ（事後）予測分布と呼びます。いま、ベイズ予測分布において $p(\theta \mid D)$ という分布を
$$ p(\theta \mid D) = \delta (\hat \theta \mid D) $$
としてみます。$\delta(a)$ とは $a$ でのみ値を持つ確率分布です。$\hat \theta$ を例えば最尤推定解に選べば、ベイズ予測分布は
$$ p(X \mid D) = p(X \mid \hat \theta) $$
と最尤推定に基づいて得られる分布と同じになります。ベイズ予測分布は $\theta$ を確率変数として捉えたことで、データに基づいて広がりを持って $\theta$ を推定し、可能性のある全ての $\theta$ を用いたアンサンブルモデルによって予測モデルを作っているということになります。
