ベイズ統計学のモデリングにおける登場人物
同時分布
ベイズ統計学で最もモデル設計者が意識しなければならないのは観測できる確率変数 $X$ と観測できない確率変数 $Z$ の同時分布です。
$$ p(X, Z) $$
ベイズ統計学でのモデリングとは主に、この同時分布にどのような仮定を与えるのかということに尽きます。
統計モデルと事前分布
例えば、 $Z$ が気温、$X$ がある計測機器を用いた気温の測定結果を表すものとしましょう。すると、実質的な気温 $Z$ の周りに、適当なばらつき $\sigma ^ 2$ で測定結果 $X$ が定まるようなことが考えられます。このモデルは、例えば下記のように平均 $Z$ 分散 $\sigma ^ 2 $ の正規分布に $X$ が従うものとして
$$ X \sim {\rm Normal} (X \mid Z, \sigma ^ 2) $$
などと書けるかもしれません。更に、 $Z$ という実際には厳密な観測が不可能な変数が、どれくらいの値を取るであろうかを確率分布 $p(Z)$ で表すこととすれば
$$ \begin{align} Z &\sim p(Z) \\ X &\sim {\rm Normal} (X \mid Z, \sigma ^ 2) \end{align} $$
なるデータの生成過程を想定することができます。通常、この $p(Z)$ のことを事前分布と呼びます。こうして、観測できるデータの統計モデルと事前分布をセットにして、今、同時分布を
$$ p(X, Z) = {\rm Normal} (X \mid Z, \sigma ^ 2)p(Z) $$
とモデル化したことになるのです（もちろん $p(Z)$ も実際には具体的な分布を仮定することになる。何も想定できないならば無情報事前分布を与える。）。 一般的に同時分布 $p(X, Z)$ に対しては
$$ p (X, Z) = p(X \mid Z)p(Z) $$
が常に成り立つため、今回は、単に統計モデル $p(X\mid Z)$ を正規分布と仮定したことがモデリングに相当することになります。
尤度関数
尤度関数は、観測できる変数 $X$ の生成過程を表す統計モデルに対して、実際に観測された値を代入して得られる関数です。尤度関数は、もはや $X$ に値を代入してしまっているので、単に未観測の $Z$ に関する関数となっています。そのことを明示するために、例えば尤度関数を $\mathcal L (Z)$ などと表記する場合もあります。
ここで、観測できる確率変数 $X$ に関しては、データが $D = \{x _ 1, ..., x _ N \}$ と集まっていることとし、尤度関数は
$$ \mathcal L (Z) = p(D \mid Z) = \prod _ i p(X = x _ i \mid Z) $$
と表記されたりします。ここで $x _ i$ は互いに（$Z$ に関して条件付きの）独立であることを想定しています。だから単にそれぞれの生成確率の積で表現できるというわけです。そうでない場合は最右辺の限りではないことに注意しましょう。
事後分布
次に重要なのが事後分布と呼ばれる モデリングした同時分布から導かれる、観測できる確率変数確定後の潜在変数の分布 です。 今、モデリングした同時分布に、尤度関数を反映してやることとすれば
$$ p(X , Z) = p(X\mid Z) p(Z) $$
に対して $X$ の実現値 $D = \{x _ 1, ..., x _ N \}$ を代入してやることにして、
$$ p(D , Z) = \prod _ i p(X = x _ i \mid Z) p(Z) $$
と書くことができます。さて、同時分布は $p(X, Z) = p(X \mid Z)p(Z) = p(Z \mid X)p(X)$ が 常に 成り立っていることに注意すれば、最左辺と最右辺に関して $p(X)$ で除してやるベイズの定理によって
$$ p(Z \mid X) = \frac{p(X,Z)}{p(X)} $$
と表記できます。この恒等式 に対して、自分がモデリングした同時分布と尤度関数を反映してやることにすれば
$$ p(Z \mid D) = \frac{p(D,Z)}{p(D)} = \frac{\prod _ i p(X = x _ i \mid Z) p(Z)}{p(D)} $$
と、観測データの実現値を得たあとの未観測変数の分布（事後分布）を書き表せるようになります。
周辺尤度
事後分布を求めることがベイズ推論における「学習」あるいは「推論」と呼ばれる処理 であるのですが、さて、上記までの手続きで、一見「学習」などと仰々しいことを言わなくとも既に事後分布は求まっているように見えます。 実際、同時分布の設計を終えた段階で、あとは数学的に成り立っている式変形をしたにすぎなく、確かに事後分布は定まっているのです。定まっているのですが、事後分布
$$ p(Z \mid D) = \frac{p(D,Z)}{p(D)} = \frac{\prod _ i p(X = x _ i \mid Z) p(Z)}{p(D)} $$
の分母 $p(D)$ を求めることの困難さによって、よほど性質の良い同時分布を使わない限り、数式で綺麗に分布を獲得することはできないのです。 この $p(D)$ のことを周辺尤度と呼びます。この周辺尤度に関して、基本に立ち返ってみてみます。
まず、確率分布で常に成り立っている関係式として、下記の式があります。
$$ p(X) = \int _ Z p(X, Z) \mathrm d Z $$
というものです。これは周辺化定理、あるいは加法定理として知られています。これは右から左に読むほうが分かりやすく、特に離散的な場合には
$$ p(X) = \sum _ Z p(X, Z) $$
と書き表されます。言っていることは単純で、2つの確率変数 $X, Z$ があるときに、同時確率 $p(X, Z)$ について $Z$ が取りうるすべてのパターンを網羅して、確率を足し合わせてしまえば、 $Z$ による確率的な効果は無くなって、$X$ だけの確率分布になります。というものです。例えば、 $Z$ は箱 $0$ と箱 $1$ のどちらを選んだかを表す確率変数で、$X$ は箱から取り出された玉の色とします。
$$ p(X, Z=0) + p(X, Z=1) $$
としてしまえば、箱 $0$ から選ばれる玉の色の確率と、箱 $1$ から選ばれる玉の色の確率を足し合わせているので、もはや玉の総数と色だけが問題になるはずで、それはすなわち $p(X)$ という確率になります。これが連続であっても常に成り立っているのですから、あえて $p(X)$ の背後に潜む潜在変数 $Z$ を仮定して、その同時分布 $p(X, Z)$ によって
$$ p(X) = \int _ Z p(X, Z) \mathrm d Z $$
と $p(X)$ を書き直すことを周辺化と呼びます。さて、周辺尤度$p(D)$は $X$ の確定値を代入して得られる確率値であるので
$$ p(D) = \int _ Z p(D, Z) \mathrm d Z = \int _ Z \prod _ i p(X = x _ i \mid Z) p(Z) \mathrm d Z $$
となっています。これは確かになかなか骨の折れそうな式になってしまっています。この積分計算が厳しいために、MCMCや変分推論などが必要になるというストーリーです。
まとめ
まとめてしまえば、ベイズ推論は同時分布
$$ p(X, Z) $$
を設計することで、ベイズの定理により
$$ p(Z \mid X) = \frac{p(X,Z)}{p(X)} $$
が具体的に定まる。分母について周辺化し
$$ p(Z \mid X) = \frac{p(X,Z)}{\int _ Z p(X, Z) \mathrm d Z} $$
と表記しておき、既に手元にある観測データ $D$ を $X$ の実現値として代入してみると
$$ p(Z \mid D) = \frac{p(D,Z)}{p(D)} = \frac{\prod _ i p(X = x _ i \mid Z) p(Z)}{\int _ Z \prod _ i p(X = x _ i \mid Z) p(Z) \mathrm d Z } $$
となって、運が良ければ数式で綺麗に求まる、多くの場合はMCMCや変分推論に頼ることになるというストーリーになります。
また、この式の重要な見方としては、 事前分布 $p(Z)$ に対して尤度関数 (観測データが反映された統計モデルを使った確率値)によって事後分布を算出しているということを、「事前分布をデータで更新した」と考えることです。 これは俗に言うベイズ更新と呼ばれるものになります（なんてことはない、単に同時分布を条件付き分布で表わして、観測データを代入しただけなのだ。そこに哲学的な話は実際のところない…）。
寄り道
最大事後確率推定（MAP）推定は、
$$ {\rm argmax} _ Z p(Z \mid D) $$
を解くことで $Z$ を点推定します。この時は厄介な周辺尤度を考える必要はないです。なぜなら、 $Z$ についての最大化問題を解く上で、周辺尤度は $Z$ を変数に持たないために影響を及ぼさないからです（積分消去されている！もはやただの定数である）。
ところでMAP推定は、そもそも $Z$ を確率変数だと考えなければ事後分布という概念が出てこないため、「ベイズ統計」を考えることで初めて出てくるアイデアのように思えてきます。その結果かは謎ですが、「MAP推定」をベイズ統計の初歩として解説する記事や本を時折目にします。それはそれで一つの考え方ではありますが、MAP推定は本質的には最尤推定に対して正則化によるペナルティを与える推定と全く変わりません。
ベイズ推論の本質はMAP推定では現れてこないと筆者は考えています。
ベイズ予測分布
ベイズ予測分布という非常に重要な分布を与えます。ベイズ予測分布とは既に手元にデータ $D$ があるとして、今後、新規のデータ $x _ {\rm new}$ が確率変数 $X$ の実現値として生ずる確率を表現する分布です。数式で表すと
$$ p(X = x _ {\rm new} \mid D) $$
と表記されます。ここで条件付き側の $D$ は表記されることもあればされないこともありますが、ベイズ予測分布の概念を考える上では非常に大事なものになります。 さて、この分布は具体的には、統計モデルと事後分布を使って下記のように構成されます。
$$ p(X = x _ {\rm new} \mid D) = \int _ Z p(X = x _ {\rm new} \mid Z)p(Z\mid D) \mathrm d Z $$
その心は、統計モデル $p(X \mid Z)$ で $x _ {new}$ の生成確率を見積もってみようとした時に、 観測できない変数 $Z$ が取りうる値を事後分布で見積もってすべて考慮した上で予測を行おうというものです。
そして、そもそも、もしも観測できる変数 $X$ に関する真の分布 $\hat p(X)$ なるものがあるとしたならば、その分布は周辺化定理によって
$$ \hat p(X) = \int _ Z \hat p(X, Z) \mathrm d Z = \int _ Z \hat p(X\mid Z) \hat p(Z) \mathrm d Z $$
となっていることが数学的に保証されています。ここで $\hat p(X \mid Z)$ と $\hat p(Z)$ はそれぞれ真の統計モデルと真の事前分布とでも言うべきでしょうか…（無理やり言葉を当てているだけだが）。そんな分布を使うことができれば $x _ {\rm new}$ の予測は
$$ \hat p(X = x _ {\rm new}) = \int _ Z \hat p(X = x _ {\rm new}, Z) \mathrm d Z = \int _ Z \hat p(X = x _ {\rm new} \mid Z) \hat p(Z) \mathrm d Z $$
となっているはずです。いずれにしてもそんなものは我々に知る由もないのですが…。 だからこそ、私達は手元にあるデータ $D$ で可能な限りベストを尽くし、事前分布 $p(Z)$ をベイズ更新して $p(Z \mid D)$ を獲得し、ドメイン知識を凝らした統計モデル $p(X \mid Z)$ を重み付けしてベイズ予測分布
$$ p(X = x _ {\rm new} \mid D) = \int _ Z p(X = x _ {\rm new} \mid Z)p(Z\mid D) \mathrm d Z $$
を構成するのです。ベイズ統計学とは、このベイズ予測分布によって真の分布を近似したり、高い予測性能を達成しようという試みだと言えます（だからMAP推定はベイズ統計で本質的ではないと言った）。
回帰の例
同時分布設計
入力 $x$ と出力 $y$ があり、その入出力関係が $y = f(x)$ と関数で表せると考えます。この関数 $f$ の形状を決めているであろうパラメータを $z$ としたものが回帰モデルです。
ここで $x$ は年齢で $y$ が年収、$z _ 1$ が年功序列に応じた昇給率で、 $z _ 0$ がボーナスであったという例を考えましょう。さて、今調査できたのは見知らぬ会社の社員の年齢とそれに紐づく年収だけであるとします。すると、$x, y$ は観測できる確率変数 $X, Y$ の実現値で、昇給率 $z _ 1$ やボーナス $z _ 0$ は潜在変数 $Z _ 0, Z _ 1$ の実現値であるということになります（ただしこちらは観測できていない）。 ここで、この問題設定では $$ y = z _ 0 + z _ 1 x $$
のような回帰モデルを考えてみることにしましょう。実際には同じ年齢であったとしてもバラツキがあるとしたならば、それも考慮してやらねばなりません。ベイズ統計学で考えるべきは
$$ p(X, Y, Z _ 0, Z _ 1) $$
です。この同時分布をどのように設計するのでしょうか。
今、昇給率 $Z _ 1$ が適当な確率分布 $p(Z _ 1)$ で仮定できるとしましょう。またボーナスも適当に $p(Z _ 0)$ で見積もっておくことにします。これらが仮に決まりさえすれば、それに対して年齢のサンプリング $p(X)$ を考慮して年収を算出しみるという手順を作ることができます。この過程は
$$ \begin{align} z _ 0, z _ 1 &\sim p(Z _ 0, Z _ 1) \\ x & \sim p(X) \\ y &\sim p(Y \mid X = x, Z _ 0 = z _ 0, Z _ 1 = z _ 1) \end{align} $$
というサンプリングをしてみたということになります。これは実は同時分布が
$$ p(X, Y, Z _ 0, Z _ 1) = p(Z _ 0, Z _ 1)p(X) p(Y \mid X, Z _ 0 , Z _ 1) $$
であると設計したということになるのです。 多くの場合、回帰モデルで実際にやりたいことは、$x$ を知っているうえで $y$ がどんな値になるかなのですから、 $p(X)$ からのサンプリングは考えずに、ある確定値を与える前提で、$X$ を条件付き側に持ってきて
$$ p(Y, Z _ 0, Z _ 1 \mid X) = \frac{p(X, Y, Z _ 0, Z _ 1)}{p(X)} = p(Z _ 0, Z _ 1) p(Y \mid X, Z _ 0 , Z _ 1) $$
という $X$ はいつでも確定値を持っているという前提での、$Y$ と潜在変数の同時分布を考えることが多いでしょう。
今、各々の確率分布がなんであるかを具体的には書いていませんが、潜在変数だけに纏わる $p(Z _ 0, Z _ 1)$ は、設計者が勝手に決めた「事前分布」と呼ばれるものです。一方で、観測変数に纏わる分布 $p(Y \mid X, Z _ 0 , Z _ 1)$ を統計モデルと呼び、通常の統計学で考慮されるのはこちらだけです（通常の統計学では観測していない値は、未知だが決定論的なパラメータとして扱う）。
当然こちらの分布も「ありえそうだ」と思う形状を設計者が勝手に決めているものです（統計モデルに実際に観測された実現値たちを定数として代入した場合に、それを尤度関数と呼ぶ。これは統計モデルを使ってデータの生起確率を計算したものであり、得られている有限のデータに対する統計モデルの当てはまり具合を計算していることになる）。
ベイズ統計モデルは、事前分布と統計モデル（尤度関数）を両方定めて設計されます。そしてそれは、実際のところ、単に「関係ありそうな変数を全て確率変数と考えた場合の同時分布」にすぎないのです。
事後分布
$p(Z _ 0, Z _ 1) $ は事前知識としてありそうな分布、あるいは何も情報を有していないということをモデルにあてがっただけのものでした。もしも、観測したデータ $D = \{x _ i, y _ i\} _ {(i = 1, \cdots, N)}$ から、この事前分布たちを修正してあげることができたら「ちゃんとデータを見て分布を考慮した」と言えます。潜在変数 $Z$ に対して、ちゃんとデータ $X, Y$ の実現値を見て考慮した条件付き分布 $p(Z \mid X, Y)$ を事後分布と呼びます。
事後分布は 常に 下記のように記せます。
$$ p(Z \mid X) = \frac{p(X, Z)}{p(X)} $$
今、観測できる確率変数は $X , Y$ の2つで、潜在変数も $Z _ 0, Z _ 1$ であったことを思い出せば
$$ p(Z _ 0, Z _ 1 \mid X, Y) = \frac{p(X,Y, Z _ 0, Z _ 1)}{p(X, Y)} $$
という形になっています（難しく考える必要はない。観測変数 $X \leftarrow \{X, Y \}$ と置き換えているだけだ。潜在変数も同様）。 そして、すでに同時分布の設計は終えているのですから、右辺の分子を 自分が設計した同時分布 に置き換えれば
$$ p(Z _ 0, Z _ 1 \mid X, Y) = \frac{p(Z _ 0, Z _ 1) p(X) p(Y \mid X, Z _ 0 , Z _ 1)}{p(X, Y)} $$
と事後分布が定まります。すなわち、統計モデルと事前分布を定めた時点で事後分布の形は数学的に決まっているのです。ところで、ここでも $X$ は常に確定値を与える前提に立てば分母に関して
$$ p(X, Y) = p(Y \mid X)p(X) $$
と変形しておけば、分母分子で $p(X)$ がキャンセルされ
$$ p(Z _ 0, Z _ 1 \mid X, Y) = \frac{p(Z _ 0, Z _ 1) p(Y \mid X, Z _ 0 , Z _ 1)}{p(Y\mid X)} $$
となります（やっていることは、自分が仮定した同時分布を代入して、あとは常に成り立っている条件付き分布の式を弄くって整理しているだけだ…！）。
さて、統計モデルに対して、実際に観測したデータ$\{x _ i, y _ i\} _ {(i = 1, \cdots, N)}$ を代入してやれば、
$$ p(Z _ 0, Z _ 1 \mid D) = \frac{p(Z _ 0, Z _ 1) \prod _ i p(X = x _ i) p(Y = y _ i \mid X = x _ i, Z _ 0 , Z _ 1) }{ \prod _ i p(Y = y _ i \mid X = x _ i)} $$
となります（ここで $D$ というのは観測したデータの具体的な値のことである。 $X, Y$ という確率変数には、もはや具体的な値が代入されている。それによって、統計モデルと事前分布によって定まるベイズモデルがデータによって具体的に計算され、事後分布の形状を更新するというわけだ。）。ともすれば、本来、人間がデータと向き合ってドメインごとに知識を凝らして行わなければならない工夫は、同時分布の設計に尽きるのです。
ベイズ予測分布
今、回帰の問題の例を続けます。 もしも $y = f(x)$ に関して、条件付き分布 $p(Y \mid X)$ を知りたいのであれば、
$$ p(Y \mid X) = \int _ {Z _ 0} \int _ {Z _ 1} p(Z _ 0, Z _ 1) p(Y \mid X, Z _ 0 , Z _ 1)\mathrm d Z _ 1 \mathrm d Z _ 0 $$
が必要なのでした。しかし、これはすでに述べたとおり叶わない願いです。そこで、まず第一の問題であった「事前分布がまともなのかどうか」という点を下記のように修正して解決します。おそらく事前分布を一発でまともなものにできるとしたら、すでにデータについて知り尽くしているか、スピノザの神かのどちらかしかありません。そういう人はベイズ統計にすがることなど無いので、通常は下記の手段を取ることになります。
$$ p(Y \mid X , D) = \int _ {Z _ 0} \int _ {Z _ 1} p(Z _ 0, Z _ 1\mid D) p(Y \mid X, Z _ 0 , Z _ 1)\mathrm d Z _ 1 \mathrm d Z _ 0 $$
ここで $D$ は予め手に入っている観測変数のデータセット $D = \{(x _ 1, y _ 1), \cdots, (x _ N, y _ N) \}$ です。このデータセット内のデータがそれぞれ $\rm i.i.d$ であると仮定すると、事後分布 $p(Z _ 0, Z _ 1 \mid D)$ は
$$ \begin{align} p(Z _ 0, Z _ 1 \mid D) &= p(Z _ 0, Z _ 1 \mid X = x _ 1, Y = y _ 1)\times \cdots \times p(Z _ 0, Z _ 1 \mid X = x _ N, Y = y _ N) \\ &= \frac{p(Z _ 0, Z _ 1) \prod _ i p(Y = y _ i \mid X = x _ i, Z _ 0 , Z _ 1)}{\int _ {Z _ 0} \int _ {Z _ 1} p(Z _ 0, Z _ 1) \prod _ i p(Y = y _ i \mid X = x _ i, Z _ 0 , Z _ 1)\mathrm d Z _ 1 \mathrm d Z _ 0} \end{align} $$
と観測変数の確定値が得られているという前提での潜在変数の条件付き分布の積で記述されることになります。この観測データすべてを使った事後分布によって周辺尤度の事前分布に相当していた部分を置き換えたものが、ベイズ予測分布であるということになります。
