毎日の天気予報を考えてみてください。明日の「気温」と「降水量」の同時分布 $p(\text{気温}, \text{降水量})$ を気象モデルが出してくれたとします。しかし、傘を持っていくかどうかを決めたいだけのとき、欲しいのは気温の情報ではなく $p(\text{降水量})$ だけです。気温という変数を「気にしないことにする」とき、内部的に何が起きているのか——これがまさに 周辺化（marginalization） が解決する問題です。
同時分布 $p(x, y)$ から $x$ だけの分布 $p(x)$ を取り出す操作、それを 周辺化 と呼び、得られた分布 $p(x)$ を 周辺分布（marginal distribution） と呼びます。離散なら和、連続なら積分で「片方の変数を消す」だけのシンプルな操作ですが、その応用範囲は驚くほど広い。
周辺化は次のような場面で必ず登場します。第一に、ベイズ推論：観測データの「証拠（evidence）」 $p(D) = \int p(D \mid \theta)\, p(\theta)\, d\theta$ はパラメータ $\theta$ を周辺化した量であり、モデル比較やマルコフ連鎖モンテカルロ法（MCMC）の正規化定数になります。第二に、隠れ変数モデル：混合ガウス分布、隠れマルコフモデル、変分オートエンコーダ（VAE）など、観測できない潜在変数 $z$ を含むモデルでは、観測の尤度 $p(x) = \sum_z p(x, z)$ を計算するのに周辺化が必要です。第三に、欠損データ：欠けた変数を周辺化してマージナル尤度に書き直すのがEMアルゴリズムの第一歩です。
本記事の内容
- 周辺化とは（一言定義） — 「興味のない変数を足して消す」操作、語源、周辺確率と周辺分布の用語整理
- 周辺化と周辺分布の直感的理解と数学的定義（離散・連続）
- 同時確率表からの計算手順と「周辺」という名前の由来
- 多変量正規分布における周辺化の閉形式解とその意味
- ベイズ推論における周辺化（証拠・周辺尤度）の役割
- 機械学習での周辺化（隠れ変数モデル、EM、VAE、ガウス過程）
- Pythonで離散・連続・ベイズ・潜在変数ケースを実装し可視化
- 周辺化が計算困難なときの近似手法（変分推論、MCMC、重点サンプリング）
この図が記事全体の地図です。左の同時分布 $p(x,y)$ から、興味のない変数（ここでは気温）をすべての値について足し合わせる（連続なら積分する）と、右の周辺分布 $p(x)$ が残ります。周辺化とは、この「足して消す」という一言に尽きる操作です。以下ではこれを離散・連続・ベイズ・機械学習の各場面で丁寧に見ていきます。
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
周辺化とは — 意味を一言で（初学者向け）
「周辺化」という言葉を初めて聞いたとき、多くの人がまず感じるのは「何を”周辺”するのか？」という戸惑いです。まずここで、最もシンプルな答えを示しておきます。
周辺化（marginalization）= 同時分布から「興味のない変数を足して消す」操作
$x$ と $y$ という2つの変数がある状況で、$y$ のことはどうでもよくて $x$ だけを見たいとき、$y$ のとりうる全ての値について確率を足し合わせる（連続なら積分する）だけで $y$ が消えて、$x$ だけの確率が残ります。
$$ p(x) = \sum_{y} p(x, y) \quad \text{（離散の場合）} $$
$$ p(x) = \int p(x, y)\, dy \quad \text{（連続の場合）} $$
たとえば「天気と気温の同時確率表があるとして、気温は関係なく”雨の確率だけ”を知りたい」——そのとき、気温の列ごとに雨の確率を全部足し合わせればよい。それが周辺化の全てです。
「周辺」という名前の由来
なぜ「周辺」と呼ぶのか、少し気になりますよね。これには具体的な由来があります。
昔から、同時確率を表（クロス集計表）で書く慣習がありました。各セルに $p(x_i, y_j)$ の値を書き並べたとき、表の右端や下端の余白（margin） に行の合計や列の合計を記入していました。その余白に現れる合計値こそが、まさに $p(x)$ や $p(y)$ なのです。
y=A y=B | 行の合計（余白）
x=1 0.20 0.10 | 0.30 ← p(x=1)
x=2 0.05 0.30 | 0.35 ← p(x=2)
x=3 0.15 0.20 | 0.35 ← p(x=3)
------+------+------+
列の合計 0.40 0.60 | 1.00
（余白）↑p(y=A) ↑p(y=B)
「表の余白（margin）に現れる確率」だから周辺確率（marginal probability）、その分布だから周辺分布（marginal distribution）。英語の “marginal” はもともと「余白の」という意味です。
図の中央の6セルが同時分布 $p(x,y)$、オレンジの右端の列が「行の合計」＝ $p(x)$、下端の行が「列の合計」＝ $p(y)$ です。まさに表の余白に周辺確率が並んでいるのが見て取れます。すべてのセルを足し上げると右下の 1.00（全確率）になることも確認できます。
周辺確率と周辺分布 — 用語の違い
似た用語が並ぶので、ここで整理しておきましょう。
| 用語 | 意味 | 例 |
|---|---|---|
| 同時分布 $p(x, y)$ | 2変数の組み合わせ全体の確率 | 「男かつA型の確率」 |
| 周辺確率 $p(x = x_0)$ | ある1点の値における周辺化した確率 | 「A型である確率 = 0.30」 |
| 周辺分布 $p(x)$ | 周辺確率を全ての $x$ の値にわたって並べたもの（分布全体） | 「血液型の分布 = [A:0.30, B:0.25, O:0.35, AB:0.10]」 |
「周辺確率」は1点の数値、「周辺分布」はその分布全体、という違いを押さえておくと以降の説明がよりクリアになります。
機械学習・ベイズでの重要性（一言サマリ）
なぜ機械学習でこれが重要かを先に一言でまとめておきます。
- ベイズ推論の「証拠」 — モデルがデータをどれだけ説明できるかを表す $p(D) = \int p(D|\theta)p(\theta)\,d\theta$ は、パラメータ $\theta$ を周辺化した量（周辺尤度）。モデル比較の核心指標です。
- 隠れ変数モデル — 混合ガウス分布やVAEでは「どのクラスタ/潜在変数から来たか」という見えない変数 $z$ を周辺化することで、観測 $x$ の尤度 $p(x) = \sum_z p(x,z)$ を計算します。
- 予測分布 — 「パラメータがこの値かも、あの値かも」という不確かさを足し合わせて（周辺化して）予測の分布 $p(y^*|x^*,D) = \int p(y^*|x^*,\theta)p(\theta|D)\,d\theta$ を作ります。
これらの詳細は後続セクションで丁寧に扱います。まず次のセクションでは、直感と数学的なイメージを固めましょう。
周辺化の直感 — 「片方の変数を気にしない」とは
周辺化を一言で表すなら「ある変数の値を全部足し合わせて、別の変数だけの確率に集約する」操作です。直感を掴むため、身近な例から始めましょう。
クラスの生徒について、$X$ を「性別（男・女）」、$Y$ を「血液型（A・B・O・AB）」という確率変数とします。クラス全体の同時分布 $p(X, Y)$ は「男でA型」「女でO型」など8通りの組み合わせそれぞれに確率を割り当てます。ところがある日、保健の先生が「血液型ごとの人数だけ知りたい」と言ってきたとします。性別という情報は要らない——というとき、私たちが内部的にやっているのは、
「血液型A型の人数 = 男のA型 + 女のA型」
という単純な足し算です。性別という変数を「消す」かわりに、その変数のすべての可能な値（男・女）について和を取る。この操作の結果として得られる $p(Y)$ が、$Y$ の 周辺分布 です。
ここで重要なのは、情報を捨てているわけではない という点です。同時分布 $p(X, Y)$ に含まれていた性別の情報は「足し合わされて」消えますが、血液型に関する正確な確率は保持されます。周辺化は「全情報のうち、興味ある変数の情報だけを抽出する」フィルターのような操作なのです。逆に言えば、一度周辺化してしまうと、消した変数についての情報は二度と元に戻せません。同時分布から周辺分布を作るのは一方通行で、$p(x, y)$ から $p(x)$ と $p(y)$ をそれぞれ得たとしても、その積 $p(x)p(y)$ は元の $p(x, y)$ に等しいとは限らない——これは独立性が成り立つときだけです。
ではなぜ「周辺」と呼ぶのか？ これは歴史的な由来があります。同時確率を表で書くと、各セルが $p(x_i, y_j)$ の値で、表の右側や下側の 余白（margin） に行の合計や列の合計を書き込む慣習がありました。その余白に書かれる行和・列和こそが、まさに $p(x)$ や $p(y)$ なのです。「表の余白に現れる確率」だから「周辺確率（marginal probability）」、その分布だから「周辺分布」。次セクションでは、この表をきちんと作りながら離散の場合の周辺化を具体的に追ってみます。
離散の周辺化 — 周辺確率表で完全に理解する
定義
$X, Y$ を離散確率変数とし、同時確率質量関数 $p(x, y) = P(X = x, Y = y)$ が与えられているとします。$X$ の周辺分布 $p(x)$ は、$Y$ のすべての値について和を取ることで得られます。
$$ \boxed{\,p(x) = \sum_y p(x, y)\,} $$
同様に、$Y$ の周辺分布は
$$ p(y) = \sum_x p(x, y) $$
です。「$Y$ を周辺化して $p(x)$ を得る」「$Y$ について和を取って $Y$ を消す」といった言い方もします。記号上は、消した変数 $y$ がシグマの中に閉じ込められて消滅する、と捉えると覚えやすいでしょう。
この式が「全確率の法則」と同等であることに気づくと、より深く理解できます。全確率の法則は、互いに排反な事象 $B_1, \dots, B_n$ で標本空間を覆うとき、
$$ P(A) = \sum_i P(A \mid B_i)\, P(B_i) = \sum_i P(A, B_i) $$
と書けます。$A = \{X = x\}$、$B_i = \{Y = y_i\}$ と置けば、周辺化の式そのものです。つまり、周辺化は全確率の法則を確率変数の言葉に翻訳したものに他なりません。
サイコロ2個の例
具体例を一つ。普通のサイコロ2個を振って、$X$ を1個目の目、$Y$ を「2個の合計」とします。$X \in \{1, \dots, 6\}$、$Y \in \{2, \dots, 12\}$ です。同時分布は
$$ p(x, y) = \begin{cases} 1/36 & y – x \in \{1, \dots, 6\} \\ 0 & \text{otherwise} \end{cases} $$
となります（2個目の目 $y – x$ が1〜6ならその組み合わせが1通り存在する）。$X$ の周辺分布を計算すると、各 $x$ に対して $y$ の取り得る値は $x+1, x+2, \dots, x+6$ の6通りなので、
$$ p(x) = \sum_{y=2}^{12} p(x, y) = 6 \times \frac{1}{36} = \frac{1}{6} $$
となり、当然ながら1個目のサイコロは一様分布 $\{1/6, 1/6, \dots, 1/6\}$ に従います。$Y$ の周辺分布のほうは、有名な「サイコロの和の分布」になります。$y = 7$ は最も組み合わせが多く $p(7) = 6/36$、$y = 2$ や $y = 12$ は端で $p(2) = p(12) = 1/36$ です。これは同時表の対角線方向に和を取った結果として自然に出てきます。
図は $Y$（2個の合計）の周辺分布 $p(y) = \sum_x p(x,y)$ です。$y=7$ で最大の $6/36$、両端の $y=2, 12$ で最小の $1/36$ という、おなじみの三角形の分布になります。1個目の目 $X$ を消して合計だけに集約した結果が、この山型として現れています。
周辺確率表（contingency table）
実データでも同じ手続きをすればよいだけです。たとえば100人にアンケートを取って、$X$ が「最も使うSNS（Twitter・Instagram・YouTube）」、$Y$ が「性別（男・女）」だったとします。クロス集計の結果が次のような周辺確率表になったとしましょう。
| 男 | 女 | 周辺 $p(x)$ | |
|---|---|---|---|
| 0.20 | 0.10 | 0.30 | |
| 0.05 | 0.30 | 0.35 | |
| YouTube | 0.15 | 0.20 | 0.35 |
| 周辺 $p(y)$ | 0.40 | 0.60 | 1.00 |
表の中央の4セル × 3行 = 6つの値が同時分布 $p(x, y)$、右端の列が $p(x) = \sum_y p(x,y)$ の周辺、下端の行が $p(y) = \sum_x p(x,y)$ の周辺です。これがまさに「表の余白に現れる確率」 = 周辺確率。すべての値を足し上げると全確率 1.0 になることも確認できます。
ここで重要な注意点があります。周辺分布だけを見ても、変数間の依存関係はわかりません。たとえば $p(\text{Instagram}) = 0.35$、$p(\text{女}) = 0.60$ という事実から「女性のうち Instagram ユーザの割合」を逆算することはできません。それを計算するには条件付き確率
$$ p(\text{Instagram} \mid \text{女}) = \frac{p(\text{Instagram}, \text{女})}{p(\text{女})} = \frac{0.30}{0.60} = 0.50 $$
が必要で、これは同時分布の値そのものを使うため、周辺化で失われた情報を取り戻す形になります。周辺化は「情報の集約」、条件付けは「情報の絞り込み」と覚えるとよいでしょう。両者は同時分布から異なる「断面」を取り出す双子の操作で、ベイズ推論ではこの両方を組み合わせて事後分布を構築します。
図のとおり、周辺化は $y$ のすべての値を足して $y$ を消す「集約」、条件付けは $y=y_0$ の断面だけを取り出して正規化する「絞り込み」です。同じ同時分布から出発しても、操作が違えば得られる分布はまったく別物になります。この2つを混同しないことが、確率モデルを正しく読むうえで重要です。
ここまでは離散の場合を見てきました。連続変数では和が積分に変わるだけですが、確率密度関数の積分が登場し、計算が一気に解析的になります。次に連続の場合を扱います。
連続の周辺化 — 積分による定式化
定義
連続確率変数 $X, Y$ の同時確率密度関数を $p(x, y)$ とすると、$X$ の周辺密度関数は
$$ \boxed{\,p(x) = \int_{-\infty}^{\infty} p(x, y)\, dy\,} $$
で定義されます。離散の総和が積分に変わっただけで、本質的な意味は同じ「$y$ のすべての可能な値について平均を取り、$y$ を消す」です。
直感的には、二変数の密度を3次元の山として描いたとき、ある $x$ の値で切った断面（$y$ について一直線に並ぶ密度のプロファイル）を集約してその断面積を求める操作、と見ることができます。$p(x)$ は「$x$ という値のところに、$y$ 軸方向にどれくらいの確率の高さが詰まっているか」を表すのです。
図は二変量正規分布を3次元の山として描き、それを左右の壁に「潰した」影が周辺分布であることを示しています。橙の曲線が $x_2$ を消した $p(x_1)$、赤の曲線が $x_1$ を消した $p(x_2)$ です。山を一方向に押しつぶす（積分する）と、その軸方向の一変量分布が壁に現れる、というのが連続の周辺化の幾何的イメージです。
二変量一様分布の例
具体例で慣れましょう。$0 \leq x, y \leq 1$ かつ $x + y \leq 1$ の三角形領域上の一様分布を考えます。三角形の面積は $1/2$ なので、密度は領域内で $p(x, y) = 2$、領域外でゼロです。$X$ の周辺密度は
$$ p(x) = \int_0^{1-x} 2\, dy = 2(1 – x), \quad 0 \leq x \leq 1 $$
となります。$x$ が大きくなるほど三角形の縦幅 $1-x$ が縮むため、周辺密度は線形に減少します。$x = 0$ で最大値2、$x = 1$ でゼロ。$\int_0^1 2(1-x)\, dx = 1$ で正規化されていることも確認できます。
多変量正規分布での周辺化
連続の世界で最も美しい結果は、多変量正規分布での周辺化が 閉形式（closed form） で求まることです。$d$ 次元正規分布 $\bm{X} = (X_1, \dots, X_d)^\top \sim \mathcal{N}(\bm{\mu}, \bm{\Sigma})$ を考え、変数を二つのブロックに分けます：
$$ \bm{X} = \begin{pmatrix} \bm{X}_1 \\ \bm{X}_2 \end{pmatrix}, \quad \bm{\mu} = \begin{pmatrix} \bm{\mu}_1 \\ \bm{\mu}_2 \end{pmatrix}, \quad \bm{\Sigma} = \begin{pmatrix} \bm{\Sigma}_{11} & \bm{\Sigma}_{12} \\ \bm{\Sigma}_{21} & \bm{\Sigma}_{22} \end{pmatrix} $$
このとき、$\bm{X}_2$ を周辺化して得られる $\bm{X}_1$ の周辺分布は次のとおりです：
$$ \boxed{\,\bm{X}_1 \sim \mathcal{N}(\bm{\mu}_1,\, \bm{\Sigma}_{11})\,} $$
驚くほどシンプルです。「$\bm{X}_1$ に対応する平均ベクトルと共分散行列の部分行列をそのまま取り出すだけ」で周辺分布が得られます。$\bm{\Sigma}_{12}$ や $\bm{\Sigma}_{22}$ は捨てられ、$\bm{X}_2$ との相関情報は消えます。これは離散の場合に表の余白で和を取ることに対応する、まさに正規分布版の「余白読み」です。
図のように、共分散行列を4ブロックに分けたとき、$\bm{X}_1$ の周辺分布は左上ブロック $\bm{\Sigma}_{11}$ をそのまま取り出すだけで得られます。積分を一切計算せず、部分行列を抜き取るだけで済むのが多変量正規の際立った性質です。捨てられる $\bm{\Sigma}_{12}, \bm{\Sigma}_{22}$ が「相関を消す」ことに対応しています。
なぜそうなるのか — 導出のスケッチ
二変量正規分布 $p(x_1, x_2)$ を $\bm{\mu} = \bm{0}$、$\bm{\Sigma} = \begin{pmatrix} \sigma_1^2 & \rho\sigma_1\sigma_2 \\ \rho\sigma_1\sigma_2 & \sigma_2^2 \end{pmatrix}$ として、$x_2$ について積分してみます。指数部を完全平方の形に書き直すと、
$$ -\frac{1}{2(1-\rho^2)}\left(\frac{x_1^2}{\sigma_1^2} – \frac{2\rho x_1 x_2}{\sigma_1 \sigma_2} + \frac{x_2^2}{\sigma_2^2}\right) $$
ここで $x_2$ について平方完成します。$x_1$ を定数とみなして $x_2$ の二次式を整理すると、
$$ \frac{1}{\sigma_2^2}\left(x_2 – \frac{\rho\sigma_2}{\sigma_1}x_1\right)^2 + \frac{(1-\rho^2) x_1^2}{\sigma_1^2} $$
の形に変形できます。第一項は $x_2$ の二次関数（中心が $x_1$ に依存）、第二項は $x_1$ だけに依存する項です。これを元の指数に戻して整理すると、
$$ p(x_1, x_2) = \underbrace{\frac{1}{\sqrt{2\pi}\sigma_1}\exp\!\left(-\frac{x_1^2}{2\sigma_1^2}\right)}_{\text{$x_1$ の周辺}} \cdot \underbrace{\frac{1}{\sqrt{2\pi(1-\rho^2)}\sigma_2}\exp\!\left(-\frac{(x_2 – \mu_{2|1})^2}{2(1-\rho^2)\sigma_2^2}\right)}_{\text{$x_2 \mid x_1$ の条件付き}} $$
と因数分解できます（$\mu_{2|1} = \frac{\rho\sigma_2}{\sigma_1}x_1$）。第二の因子は $x_2$ について正規分布の密度なので、積分は1。残るのは第一の因子で、これがまさに $\mathcal{N}(0, \sigma_1^2)$ の密度です。
この導出が一般の多変量に拡張されると、「平均と共分散行列の部分行列をそのまま使う」という結論に到達します。詳しくは 多変量正規分布の周辺化 で扱います。
連続の場合は積分が解析的に解ければ閉形式で求まりますが、解けない場合は数値積分や近似が必要になります。次に、応用上もっとも重要な「ベイズ推論における周辺化」を見ていきましょう。
ベイズ推論における周辺化 — 証拠と周辺尤度
ベイズの定理に潜む周辺化
ベイズの定理を書き下すと
$$ p(\theta \mid D) = \frac{p(D \mid \theta)\, p(\theta)}{p(D)} $$
ここで $\theta$ がパラメータ、$D$ が観測データです。分母の $p(D)$ は 証拠（evidence） あるいは 周辺尤度（marginal likelihood） と呼ばれ、
$$ p(D) = \int p(D, \theta)\, d\theta = \int p(D \mid \theta)\, p(\theta)\, d\theta $$
と書けます。これはまさに $\theta$ について周辺化した量です。「データ $D$ がどれくらいモデル全体から見て自然か」を、すべての可能なパラメータについて尤度を平均したもの。「証拠」という名前は、「観測データがそのモデルを支持する強さの証拠」という意味です。
モデル比較
二つのモデル $M_1, M_2$ を比較したいとき、ベイズファクター
$$ \mathrm{BF}_{12} = \frac{p(D \mid M_1)}{p(D \mid M_2)} $$
は両モデルの周辺尤度の比です。$p(D \mid M_k) = \int p(D \mid \theta, M_k)\, p(\theta \mid M_k)\, d\theta$ は、モデル $k$ のパラメータ事前分布で平均したデータ尤度。周辺尤度はモデル比較の根本指標 なのです。最大尤度（MLE）が「最良パラメータでの当てはまり」を測るのに対し、周辺尤度は「全パラメータを考慮したモデル全体の予測能力」を測ります。複雑すぎるモデルは事前分布が広がりすぎて周辺尤度が下がる（オッカムの剃刀の自動実装）、というのもこの周辺化のおかげです。
予測分布
新しい入力 $x^*$ に対する予測 $y^*$ を作るときも、パラメータの不確かさを周辺化します：
$$ p(y^* \mid x^*, D) = \int p(y^* \mid x^*, \theta)\, p(\theta \mid D)\, d\theta $$
これは「事後分布 $p(\theta \mid D)$ に従って $\theta$ をサンプルし、各 $\theta$ での予測 $p(y^* \mid x^*, \theta)$ の平均を取る」操作です。点推定（MLEやMAPで一つの $\theta$ に固定）と異なり、パラメータの不確かさが予測の不確かさに正しく反映されます。
共役事前分布での閉形式
事前分布と尤度が共役関係にあれば、周辺尤度は閉形式で求まります。たとえば、ベルヌーイ試行 $x_i \sim \mathrm{Bernoulli}(\theta)$ に対してベータ事前分布 $\theta \sim \mathrm{Beta}(\alpha, \beta)$ を置いたとき、$n$ 回中 $k$ 回成功したデータの周辺尤度は
$$ p(D) = \int_0^1 \binom{n}{k}\theta^k(1-\theta)^{n-k} \cdot \frac{\theta^{\alpha-1}(1-\theta)^{\beta-1}}{B(\alpha,\beta)}\, d\theta = \binom{n}{k}\frac{B(\alpha+k, \beta+n-k)}{B(\alpha,\beta)} $$
と、ベータ関数 $B(\cdot, \cdot)$ を使って簡潔に表せます。詳しくは ベイズ推定とは？仕組みについてざっくり解説 を参照してください。
なぜ周辺尤度の計算は難しいのか
理論的には美しい周辺尤度ですが、実用上は計算が非常に困難なことが多いです。$\theta$ が高次元（深層学習なら数百万次元）になると、積分 $\int p(D \mid \theta)\, p(\theta)\, d\theta$ は数値的に手に負えません。共役性が成り立たない一般のモデルでは、閉形式が存在せず、変分推論やMCMC（マルコフ連鎖モンテカルロ）といった近似手法が不可欠になります。
ベイズ推論での周辺化は「パラメータの不確かさを統合する」操作でしたが、機械学習にはもう一種類の重要な周辺化があります——観測できない隠れ変数の周辺化 です。
機械学習での周辺化 — 隠れ変数とEMアルゴリズム
隠れ変数モデル
機械学習の多くのモデルは、観測変数 $x$ と 隠れ変数（潜在変数、latent variable） $z$ を持ちます。$z$ は「データが実際にはどのクラスタから来たか」「文書のトピックは何か」「画像の本質的な特徴は何か」といった、観測できないが構造を生成する内部状態です。同時分布 $p(x, z \mid \theta)$ を直接モデル化するのが普通ですが、対数尤度を計算するには $z$ を周辺化する必要があります：
$$ p(x \mid \theta) = \sum_z p(x, z \mid \theta) \quad \text{or} \quad \int p(x, z \mid \theta)\, dz $$
これがいわゆる 観測尤度（marginal likelihood） で、最尤推定で最大化したい量です。
混合ガウス分布
代表例として混合ガウス分布（GMM）を見ましょう。$K$ 個の正規分布の混合として
$$ p(x \mid \theta) = \sum_{k=1}^K \pi_k\, \mathcal{N}(x \mid \mu_k, \Sigma_k) $$
と書けます。これは隠れ変数 $z \in \{1, \dots, K\}$（どのクラスタから来たか）を導入し、$p(z = k) = \pi_k$、$p(x \mid z = k) = \mathcal{N}(x \mid \mu_k, \Sigma_k)$ としたときの同時分布 $p(x, z) = p(z) p(x \mid z)$ を $z$ について周辺化したものです。
$$ p(x) = \sum_{k=1}^K p(x, z=k) = \sum_{k=1}^K \pi_k\, \mathcal{N}(x \mid \mu_k, \Sigma_k) $$
つまり「データが属するクラスタという情報を消す」操作。直感的には、各正規分布の山を重み付きで合成した分布。$K$ が小さければ和は数項なので扱いやすいですが、隠れマルコフモデル（HMM）のように $z$ が系列なら、状態空間のサイズは指数的に爆発します。
EMアルゴリズム
直接 $\log p(x \mid \theta)$ を最大化するのは、和の中に対数があるため難しい。そこで EMアルゴリズム は周辺尤度の下界（変分下界、ELBO）を最大化する迂回路を取ります：
$$ \log p(x \mid \theta) \geq \mathbb{E}_{q(z)}[\log p(x, z \mid \theta)] – \mathbb{E}_{q(z)}[\log q(z)] $$
$q(z)$ は $z$ の仮想的な分布。Eステップ で $q(z) = p(z \mid x, \theta)$ に更新（下界がタイトになる）、Mステップ で $\theta$ を更新（下界を最大化）。この反復で対数周辺尤度が単調に増加します。周辺化を直接やる代わりに、$q(z)$ を介して間接的に隠れ変数を扱う技法です。
変分オートエンコーダ（VAE）
深層学習における周辺化の象徴的応用がVAEです。画像 $x$ の生成過程を、潜在変数 $z$（低次元の意味的特徴）からニューラルネットでデコードするモデル $p_\theta(x, z) = p(z) p_\theta(x \mid z)$ として書きます。観測尤度
$$ p_\theta(x) = \int p(z) p_\theta(x \mid z)\, dz $$
は高次元の積分で計算不可能。そこで近似事後 $q_\phi(z \mid x)$（エンコーダ）を学習し、ELBOを最大化します：
$$ \mathcal{L}(\theta, \phi; x) = \mathbb{E}_{q_\phi(z|x)}[\log p_\theta(x \mid z)] – \mathrm{KL}(q_\phi(z \mid x) \| p(z)) $$
ELBOの最大化は周辺尤度の下界の最大化、すなわち間接的な周辺化操作です。詳細は KL距離・JSダイバージェンスとは を参照してください。
ガウス過程回帰
ガウス過程（GP）は無限次元のパラメータ（関数空間）を持つベイズモデルですが、多変量正規分布の周辺化の性質のおかげで実用化されています。観測点 $\bm{x}$ と新規点 $\bm{x}^*$ の関数値 $f$ の同時分布が多変量正規分布
$$ \begin{pmatrix} \bm{f} \\ \bm{f}^* \end{pmatrix} \sim \mathcal{N}\!\left(\bm{0}, \begin{pmatrix} \bm{K} & \bm{K}_* \\ \bm{K}_*^\top & \bm{K}_{**} \end{pmatrix}\right) $$
となるので、新規点 $\bm{f}^*$ について周辺化した予測分布は、再び正規分布の閉形式で得られます。無限次元の積分を計算する代わりに、有限次元の正規分布の部分行列を取り出すだけ——これが多変量正規の周辺化の素晴らしさです。
機械学習における周辺化の重要性は十分わかったと思います。次に、これらを実際にPythonで動かして確認してみましょう。
Python実装 — 周辺化を多面的に確認する
離散の周辺確率表
まずは離散の周辺化を実装します。同時分布 $p(x, y)$ を2次元 numpy 配列で表現し、np.sum の axis 引数で和を取るだけです。
import numpy as np
import matplotlib.pyplot as plt
# 同時分布 p(x, y) を行列で表現
# 行: SNS (Twitter, Instagram, YouTube), 列: 性別 (男, 女)
joint = np.array([
[0.20, 0.10], # Twitter
[0.05, 0.30], # Instagram
[0.15, 0.20], # YouTube
])
sns_labels = ['Twitter', 'Instagram', 'YouTube']
gender_labels = ['Male', 'Female']
# 周辺化: 性別について和を取る → p(SNS)
p_x = joint.sum(axis=1) # axis=1 で列方向に和
# 周辺化: SNSについて和を取る → p(性別)
p_y = joint.sum(axis=0) # axis=0 で行方向に和
print("同時分布 p(x, y):")
print(joint)
print(f"\n周辺分布 p(SNS) = {p_x} (合計 {p_x.sum():.2f})")
print(f"周辺分布 p(性別) = {p_y} (合計 {p_y.sum():.2f})")
print(f"全確率: {joint.sum():.2f}")
# 周辺確率表として可視化
fig, ax = plt.subplots(figsize=(7, 5))
im = ax.imshow(joint, cmap='Blues', vmin=0, vmax=0.35)
for i in range(joint.shape[0]):
for j in range(joint.shape[1]):
ax.text(j, i, f'{joint[i, j]:.2f}', ha='center', va='center', fontsize=12)
ax.set_xticks(range(2))
ax.set_xticklabels([f'{g}\np(y)={p:.2f}' for g, p in zip(gender_labels, p_y)])
ax.set_yticks(range(3))
ax.set_yticklabels([f'{s}\np(x)={p:.2f}' for s, p in zip(sns_labels, p_x)])
ax.set_title('Joint distribution p(x, y) with marginals')
plt.colorbar(im, ax=ax, label='Probability')
plt.tight_layout()
plt.savefig('marginalization_table.png', dpi=150, bbox_inches='tight')
plt.show()
このコードを実行すると、p(SNS) = [0.30, 0.35, 0.35]、p(性別) = [0.40, 0.60] が出力されます。同時分布の各セルが周辺分布の対応する行・列に集約されていることが確認できます。可視化されたヒートマップでは、軸ラベルに周辺確率の値が併記され、「表の余白に現れる確率」という名前の由来が直接見て取れます。同時分布の値をすべて足し上げると 1.00 で、確率の正規化条件 $\sum_{x,y} p(x,y) = 1$ が成立しています。
二変量正規分布の周辺化
連続の場合に移ります。二変量正規分布の周辺化が部分行列の取り出しになることを、サンプリングと密度プロットで確認します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal, norm
# 二変量正規分布のパラメータ
mu = np.array([1.0, 2.0])
sigma1, sigma2 = 1.5, 1.0
rho = 0.7 # 相関係数
Sigma = np.array([
[sigma1**2, rho * sigma1 * sigma2],
[rho * sigma1 * sigma2, sigma2**2]
])
print(f"平均 mu = {mu}")
print(f"共分散行列 Sigma =\n{Sigma}")
# 同時分布 p(x1, x2)
x1 = np.linspace(-3, 5, 200)
x2 = np.linspace(-2, 6, 200)
X1, X2 = np.meshgrid(x1, x2)
pos = np.dstack([X1, X2])
joint_pdf = multivariate_normal(mu, Sigma).pdf(pos)
# x2 を周辺化 → x1 の周辺分布は N(mu[0], sigma1^2)
# 数値積分でも確認
p_x1_numeric = joint_pdf.sum(axis=0) * (x2[1] - x2[0])
p_x1_theory = norm(mu[0], sigma1).pdf(x1)
# x1 を周辺化 → x2 の周辺分布は N(mu[1], sigma2^2)
p_x2_numeric = joint_pdf.sum(axis=1) * (x1[1] - x1[0])
p_x2_theory = norm(mu[1], sigma2).pdf(x2)
# 可視化: 同時分布 + 両軸に周辺分布
fig = plt.figure(figsize=(10, 8))
gs = fig.add_gridspec(2, 2, width_ratios=[3, 1], height_ratios=[1, 3],
hspace=0.05, wspace=0.05)
ax_joint = fig.add_subplot(gs[1, 0])
ax_marg_x = fig.add_subplot(gs[0, 0], sharex=ax_joint)
ax_marg_y = fig.add_subplot(gs[1, 1], sharey=ax_joint)
# 同時分布
ax_joint.contourf(X1, X2, joint_pdf, levels=15, cmap='viridis')
ax_joint.set_xlabel('$x_1$')
ax_joint.set_ylabel('$x_2$')
# x1 の周辺分布(上)
ax_marg_x.plot(x1, p_x1_numeric, 'b-', lw=2, label='Numerical $\\int p(x_1,x_2)dx_2$')
ax_marg_x.plot(x1, p_x1_theory, 'r--', lw=1.5, label='Theory $\\mathcal{N}(\\mu_1,\\sigma_1^2)$')
ax_marg_x.legend(fontsize=8)
ax_marg_x.set_title('Marginal of $x_1$')
plt.setp(ax_marg_x.get_xticklabels(), visible=False)
# x2 の周辺分布(右)
ax_marg_y.plot(p_x2_numeric, x2, 'b-', lw=2)
ax_marg_y.plot(p_x2_theory, x2, 'r--', lw=1.5)
ax_marg_y.set_title('Marginal\nof $x_2$', fontsize=10)
plt.setp(ax_marg_y.get_yticklabels(), visible=False)
plt.savefig('bivariate_marginalization.png', dpi=150, bbox_inches='tight')
plt.show()
# 最大誤差を確認
err_x1 = np.max(np.abs(p_x1_numeric - p_x1_theory))
err_x2 = np.max(np.abs(p_x2_numeric - p_x2_theory))
print(f"x1 周辺分布の最大誤差(数値 vs 理論): {err_x1:.2e}")
print(f"x2 周辺分布の最大誤差(数値 vs 理論): {err_x2:.2e}")
このグラフから、二変量正規分布の周辺化の本質がはっきり読み取れます。中央の等高線は同時分布で、相関 $\rho = 0.7$ により楕円が右上がりに傾いています。上下の余白に描かれた周辺分布（数値積分結果と理論曲線）は完全に一致し、$x_2$ を消した結果が単純な一変量正規分布 $\mathcal{N}(1.0, 1.5^2)$ になることを確認できます。理論で示した「部分行列を取り出すだけ」という結論が、数値的にも誤差 $10^{-4}$ 以下で成立しています。相関情報 $\rho$ は周辺化の過程で消えており、軸方向の幅 $\sigma_1, \sigma_2$ だけが残っていることに注目してください。
中央の等高線が同時分布（相関 $\rho=0.7$ で楕円が傾く）、上と右の余白が周辺分布です。数値積分（シアン）と理論曲線（破線）がぴったり重なり、$x_2$ を消した結果が単純な一変量正規 $\mathcal{N}(1.0, 1.5^2)$ になることが確認できます。
ベイズ推論の周辺尤度
ベイズの定理の分母である周辺尤度を、ベルヌーイ-ベータモデルで実装してみます。
import numpy as np
from scipy.special import beta as beta_func
from scipy.special import comb
import matplotlib.pyplot as plt
def marginal_likelihood_beta_bernoulli(k, n, alpha, beta):
"""ベルヌーイ-ベータの周辺尤度 p(D) = C(n,k) * B(α+k, β+n-k) / B(α, β)"""
return comb(n, k) * beta_func(alpha + k, beta + n - k) / beta_func(alpha, beta)
def marginal_likelihood_numerical(k, n, alpha, beta, num_theta=1000):
"""数値積分で周辺尤度を計算（検証用）"""
theta = np.linspace(1e-6, 1 - 1e-6, num_theta)
likelihood = comb(n, k) * theta**k * (1 - theta)**(n - k)
prior = theta**(alpha - 1) * (1 - theta)**(beta - 1) / beta_func(alpha, beta)
return np.trapz(likelihood * prior, theta)
# データ: 10回中7回成功
k, n = 7, 10
# 異なる事前分布で周辺尤度を比較
priors = [(1, 1), (5, 5), (1, 5), (10, 2)]
prior_names = ['Uniform Beta(1,1)', 'Symmetric Beta(5,5)',
'Skeptical Beta(1,5)', 'Confident Beta(10,2)']
print(f"データ: {n}回中{k}回成功\n")
print("事前分布 | 周辺尤度(閉形式) | 周辺尤度(数値積分)")
print("-" * 65)
for (a, b), name in zip(priors, prior_names):
p_D_closed = marginal_likelihood_beta_bernoulli(k, n, a, b)
p_D_numeric = marginal_likelihood_numerical(k, n, a, b)
print(f"{name:24s}| {p_D_closed:.6f} | {p_D_numeric:.6f}")
# 事前分布と尤度と事後分布を可視化
theta = np.linspace(0, 1, 200)
fig, axes = plt.subplots(1, 4, figsize=(16, 4), sharey=True)
for ax, (a, b), name in zip(axes, priors, prior_names):
likelihood = comb(n, k) * theta**k * (1 - theta)**(n - k)
prior = theta**(a - 1) * (1 - theta)**(b - 1) / beta_func(a, b)
p_D = marginal_likelihood_beta_bernoulli(k, n, a, b)
posterior = likelihood * prior / p_D
ax.plot(theta, prior, 'g-', label='Prior')
ax.plot(theta, likelihood / likelihood.max() * prior.max(),
'b--', label='Likelihood (scaled)')
ax.plot(theta, posterior, 'r-', lw=2, label='Posterior')
ax.fill_between(theta, 0, likelihood * prior, alpha=0.2, color='gray',
label=f'p(D) = {p_D:.4f}')
ax.set_title(name, fontsize=10)
ax.set_xlabel('$\\theta$')
ax.legend(fontsize=7)
ax.grid(True, alpha=0.3)
axes[0].set_ylabel('Density')
plt.suptitle('Marginal likelihood under different priors (k=7, n=10)', y=1.02)
plt.tight_layout()
plt.savefig('marginal_likelihood.png', dpi=150, bbox_inches='tight')
plt.show()
このプロットから、事前分布の違いが周辺尤度（証拠）にどう影響するかが読み取れます。事前分布が観測データと整合する場合（たとえば Beta(10,2) は「成功率が高い」と事前に予想しており、k=7/n=10 と相性が良い）、周辺尤度は大きくなります。逆に、観測と矛盾する事前（Beta(1,5) は失敗が多いと予想）では、周辺尤度は小さくなります。閉形式と数値積分の結果は $10^{-5}$ 以下の精度で一致しており、解析解の正しさが確認できます。この値はベイズファクターによるモデル比較で「どの事前分布が観測を最もよく説明するか」を定量化する基礎になります。
図は4つの事前分布での事前（緑）・尤度（青破線）・事後（赤）と、灰色の面積で表した周辺尤度 $p(D)$ を並べたものです。観測（k=7/n=10）と整合する確信的な事前 Beta(10,2) では証拠が大きく、矛盾する懐疑的な事前 Beta(1,5) では小さくなります。周辺尤度が「データがそのモデルをどれだけ支持するか」を測っていることが見て取れます。
隠れ変数の周辺化（混合ガウス）
最後に、混合ガウス分布での隠れ変数 $z$ の周辺化を実装します。各成分の重み付き和が観測分布になることを確かめます。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
# 混合ガウスのパラメータ: K=3成分
pi = np.array([0.3, 0.5, 0.2]) # 混合重み
mus = np.array([-2.0, 1.0, 4.0]) # 各成分の平均
sigmas = np.array([0.6, 0.8, 1.2]) # 各成分の標準偏差
# 観測尤度 p(x) = sum_k pi_k * N(x | mu_k, sigma_k^2)
x = np.linspace(-5, 8, 500)
components = np.array([pi[k] * norm(mus[k], sigmas[k]).pdf(x) for k in range(3)])
p_x = components.sum(axis=0) # z を周辺化 = 成分について和
# サンプリングで検証: 各サンプルで z をサンプル → x をサンプル
np.random.seed(42)
n_samples = 10000
z_samples = np.random.choice(3, size=n_samples, p=pi)
x_samples = np.array([np.random.normal(mus[z], sigmas[z]) for z in z_samples])
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
# 左: 同時分布 p(x, z) を視覚化 + 周辺化の様子
colors = ['tab:blue', 'tab:orange', 'tab:green']
for k in range(3):
axes[0].fill_between(x, components[k], alpha=0.4, color=colors[k],
label=f'$\\pi_{k}\\mathcal{{N}}(\\mu_{k},\\sigma_{k}^2)$')
axes[0].plot(x, p_x, 'k-', lw=2.5, label='$p(x) = \\sum_z p(x, z)$')
axes[0].set_xlabel('$x$')
axes[0].set_ylabel('Density')
axes[0].set_title('Marginal $p(x)$ as sum over hidden $z$')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 右: サンプリングと理論密度の比較
axes[1].hist(x_samples, bins=60, density=True, alpha=0.5, color='gray',
label=f'{n_samples} samples')
axes[1].plot(x, p_x, 'r-', lw=2.5, label='Theory $p(x)$')
axes[1].set_xlabel('$x$')
axes[1].set_ylabel('Density')
axes[1].set_title('Sampling: $z\\sim\\mathrm{Cat}(\\pi)$, $x\\sim\\mathcal{N}(\\mu_z,\\sigma_z^2)$')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('mixture_marginalization.png', dpi=150, bbox_inches='tight')
plt.show()
print(f"理論密度の積分: {np.trapz(p_x, x):.4f} (should be 1.0)")
print(f"サンプル平均: {x_samples.mean():.3f}, 理論平均: {(pi*mus).sum():.3f}")
左のプロットは「隠れ変数を周辺化する」操作を視覚化したものです。3つのガウス成分（色付き）が個別に存在し、それらを重み付きで合成（和を取る）した黒線が観測分布 $p(x) = \sum_z p(x, z)$ になります。$z$ という情報を捨てて $x$ だけの分布に集約する、まさに周辺化の操作です。右のプロットでは、$z$ をカテゴリカル分布からサンプリングして $x$ を生成する手続き（サンプリングによる周辺化）の結果が、理論密度と一致することが確認できます。サンプル平均は理論値 $\sum_k \pi_k \mu_k = 0.3 \cdot (-2) + 0.5 \cdot 1 + 0.2 \cdot 4 = 0.7$ に近づき、十分なサンプル数で周辺化が再現できることが見て取れます。
図の色付きの3つの山が各ガウス成分 $\pi_k \mathcal{N}(\mu_k, \sigma_k^2)$、黒の太線がそれらを足し合わせた観測分布 $p(x) = \sum_z p(x,z)$ です。「どの山から来たか」という隠れ変数 $z$ を消して $x$ だけの分布に集約する——これが隠れ変数モデルにおける周辺化そのものです。
計算困難な周辺化の近似 — 変分推論・MCMC・重点サンプリング
ここまで見てきたとおり、周辺化は「和か積分」のシンプルな操作です。しかし実用では、その計算自体が極端に難しい状況が頻繁に発生します。
困難になる典型例
- 高次元積分: $\theta$ が数百〜数百万次元（深層モデルのパラメータ）の場合、グリッド積分は次元の呪いで破綻します。
- 解析的に解けない密度: 事前分布と尤度が共役でないと閉形式が存在せず、被積分関数が手に負えない形になる。
- 離散変数の組み合わせ爆発: HMMの $T$ 時刻の隠れ状態は $K^T$ 通り、ベイズネットの全配置は指数爆発します。
これらに対する近似手法を概観します。
図は、計算不可能な周辺化積分に対する代表的な近似手法を整理したものです。サンプル平均で攻めるモンテカルロ積分・MCMC、最適化に帰着させる変分推論（VI）、モード周りでガウス近似するLaplace近似——問題の次元・共役性・滑らかさに応じて使い分けます。「周辺化」という一つの操作の背後に、現代統計・機械学習の実用技法の大半が連なっています。
モンテカルロ積分
最も汎用的なのが モンテカルロ積分と重点サンプリング。$p(\theta)$ から $\theta_i$ を $M$ 個サンプリングし、
$$ \int f(\theta)\, p(\theta)\, d\theta \approx \frac{1}{M}\sum_{i=1}^M f(\theta_i) $$
と平均で近似します。誤差は $O(1/\sqrt{M})$ で次元によらないため、高次元でも使えます。ただし $p(\theta)$ から効率的にサンプリングできない場合は、別の提案分布 $q(\theta)$ から取り、
$$ \int f(\theta) p(\theta)\, d\theta \approx \frac{1}{M}\sum_i f(\theta_i)\frac{p(\theta_i)}{q(\theta_i)} $$
と重みを補正する 重点サンプリング を使います。
MCMC（マルコフ連鎖モンテカルロ）
事後分布 $p(\theta \mid D)$ から直接サンプリングすることで、$\int g(\theta) p(\theta \mid D) d\theta \approx \frac{1}{M}\sum_i g(\theta_i)$ で予測や期待値を推定する手法。Metropolis-Hastings、Gibbsサンプリング、Hamiltonian Monte Carlo（HMC）、NUTSなど多彩なアルゴリズムが開発されています。周辺尤度 $p(D)$ そのものは普通のMCMCでは直接求まらない ことに注意（事後分布の正規化定数だから）。周辺尤度を欲しい場合は、ブリッジサンプリング・パスサンプリング・ネストサンプリングといった特殊技法が必要です。
変分推論（VI）
事後分布 $p(\theta \mid D)$ を解析的に扱いやすい分布族 $q_\phi(\theta)$ で近似し、KLダイバージェンス $\mathrm{KL}(q_\phi \| p)$ を最小化することで周辺化の代わりにする手法。
$$ \log p(D) = \mathcal{L}(\phi) + \mathrm{KL}(q_\phi(\theta) \| p(\theta \mid D)) $$
ELBO $\mathcal{L}(\phi) = \mathbb{E}_{q_\phi}[\log p(D, \theta) – \log q_\phi(\theta)]$ を最大化すると、KLが最小化され同時に対数周辺尤度の下界が得られます。MCMCより高速（最適化問題に帰着）ですが、$q_\phi$ の選び方によりバイアスが残ります。VAEはこのVIをニューラルネットで実装した特殊形です。
Laplace近似
事後分布のモード $\hat\theta_{\mathrm{MAP}}$ の周りで対数事後分布を二次近似し、事後を正規分布で代用する古典的手法。
$$ p(D) \approx p(D \mid \hat\theta)\, p(\hat\theta)\, \frac{(2\pi)^{d/2}}{|\bm{H}|^{1/2}} $$
ここで $\bm{H}$ はヘシアン行列。簡単で計算量が少ない反面、事後が単峰でガウシアン近似が良い場合のみ有効です。BIC（ベイズ情報量規準）はこのLaplace近似の漸近的形式から導かれます。
実用では、問題の構造（次元数、滑らかさ、共役性）に応じてこれらを使い分けます。深層ベイズなら確率的変分推論（SVI）、複雑な階層モデルならNUTS（PyMC、Stan）、リアルタイム性能が必要ならLaplace近似——といった具合です。「周辺化」という一つの数学的操作の背後に、現代統計・機械学習の実用技法の大半が連なっていることがわかります。
まとめ
本記事では、周辺化（marginalization）と周辺分布について、定義・直感・応用・実装の各側面から解説しました。
- 基本操作: 同時分布 $p(x, y)$ から $y$ について和（離散 $\sum_y p(x,y)$）または積分（連続 $\int p(x,y)\,dy$）を取って、$y$ の情報を消し $p(x)$ を得る。「表の余白に現れる確率」だから「周辺確率」。
- 離散の例: 周辺確率表で行和・列和を取るだけ。サイコロ和の分布もこの操作で導出される。情報を捨てるが、興味ある変数の確率は完全に保存される。
- 連続の例: 多変量正規分布では、平均ベクトルと共分散行列の対応部分行列をそのまま取り出すだけで周辺分布が得られる（閉形式）。
- ベイズ推論: 周辺尤度 $p(D) = \int p(D \mid \theta) p(\theta)\, d\theta$ はパラメータ $\theta$ を周辺化した量で、モデル比較（ベイズファクター）と予測分布の基礎。
- 機械学習: 隠れ変数モデル（GMM、HMM、VAE）では潜在 $z$ の周辺化が観測尤度を生む。EMアルゴリズム・ELBO最大化・ガウス過程の予測など、現代機械学習の核心技法はすべて周辺化を扱う技術。
- 計算が難しい場合: モンテカルロ積分、MCMC、変分推論、Laplace近似で近似する。高次元・非共役・組み合わせ爆発などが困難の原因。
周辺化を理解すると、確率モデルの「どの変数を保持し、どの変数を消すか」という設計思想が見えてきます。ベイズ的予測・潜在変数モデル・モデル選択など、応用先は驚くほど広く、確率論を実用に橋渡しする中核操作と言えます。
次のステップとして、以下の記事も参考にしてください。