第7章確率モデルと生成モデル（機械学習の数理）

この章の目標

エントロピー・交差エントロピー・KL ダイバージェンス・相互情報量の基本性質とギブスの不等式を証明し、最尤推定が経験分布との KL ダイバージェンスの最小化であることを示せる

混合ガウスモデルの EM アルゴリズムを導いて対数尤度の単調性を証明し、局所解と尤度の非有界性という限界を説明できる

$\log p(x) = \operatorname{ELBO}(q) + \operatorname{KL}(q \Vert p(\cdot \mid x))$ を導き、平均場近似の更新式と、EM が変分推論の特別な場合であることを説明できる

隠れマルコフモデルの前向きアルゴリズムを導き、計算量が $O(TK^2)$ であることを示せる

変分オートエンコーダ・拡散モデル・大規模言語モデルを、確率モデルと学習の目的関数の言葉で読める

前提：第1章、22-statistics 第1章（条件付き分布・多変量正規分布）、22-statistics 第3章（最尤推定）。7.5 節では 11-probability 第6章のマルコフ連鎖を、7.6 節では第2章のソフトマックス回帰と第5章のニューラルネットワークを使う。22-statistics 第7章（事後分布・マルコフ連鎖モンテカルロ法）と比べるとよい。

工場のセンサーの値がふだんの分布から外れたら異常を疑う。購入額の分布に山が 2 つあれば、性質の違う客が混ざっているのだろう。設備の劣化は直接見えず、アラームの記録から推し量るしかない。いずれも、データ $x$ の分布 $p(x)$ を、観測できない潜在変数 (latent variable) も含めてモデル化する問題である。データの生まれ方を確率分布で記述し、新しいデータを生成できるモデルを生成モデル (generative model) という。

本章では KL ダイバージェンスと最尤推定（7.1・7.2 節）、EM アルゴリズム（7.3 節）、変分推論（7.4 節）、隠れマルコフモデル（7.5 節）を学び、近年の生成モデルを同じ言葉で読む（7.6 節）。 $\log$ は自然対数、「分布 $p$ 」は確率関数または密度を指し、和は連続の場合は積分に読み替える。測度論的な細部は 11-probability に譲る。

7.1 エントロピーと KL ダイバージェンス

定義 7.1（エントロピー・交差エントロピー・KL ダイバージェンス） $p, q$ を有限集合または可算集合上の確率分布とし、 $A = \lbrace x \mid p(x) > 0 \rbrace$ とする。
$H(p) = -\sum_{x \in A} p(x)\log p(x), \qquad H(p, q) = -\sum_{x \in A} p(x)\log q(x), \qquad \operatorname{KL}(p \Vert q) = \sum_{x \in A} p(x)\log\frac{p(x)}{q(x)}$
をそれぞれエントロピー (entropy)、交差エントロピー (cross-entropy)、KL ダイバージェンス（カルバック–ライブラー情報量, Kullback–Leibler divergence）という。 $q(x) = 0$ となる $x \in A$ があれば後の 2 つは $\infty$ とする。密度でも和を $A$ 上の積分に変えて同様に定める（ $A$ 上で $q = 0$ となる部分の確率が正なら後の 2 つは $\infty$ ）。 $X \sim p$ のとき $H(X) = H(p)$ とも書く。

$-\log p(x)$ は起こりにくい値ほど大きい「驚き」で、エントロピーはその平均（分布の不確かさ）である。KL ダイバージェンスの項は負にもなるが、下の $\log u \leq u - 1$ より $p\log\frac{p}{q} \geq p - q \geq -q$ なので、負の項の絶対値の和は $1$ 以下であり、和は $(-\infty, \infty]$ の値として定まる。 $H(p) < \infty$ なら $H(p, q) = H(p) + \operatorname{KL}(p \Vert q)$ である。密度のエントロピー（微分エントロピー）は負にもなる。

定理 7.2（ギブスの不等式, Gibbs' inequality） $\operatorname{KL}(p \Vert q) \geq 0$ であり、等号は $p = q$ のときに限る。密度の場合も $\operatorname{KL}(p \Vert q) \geq 0$ で、等号は $p$ と $q$ がほとんど至るところ（ルベーグ測度 $0$ の集合を除いて）等しいときに限る。

証明. $g(u) = u - 1 - \log u$ （ $u > 0$ ）は $g'(u) = 1 - 1/u$ が $u = 1$ の前後で負から正に変わるので、最小値 $g(1) = 0$ をとる。すなわち $\log u \leq u - 1$ で、等号は $u = 1$ に限る。 $q = 0$ となる $A$ の点があれば $\operatorname{KL} = \infty$ なので、 $A$ 上 $q > 0$ とする。 $p\log\frac{q}{p} = (q - p) - p \cdot g(q/p)$ より

-\operatorname{KL}(p \Vert q) = \sum_{x \in A}\bigl(q(x) - p(x)\bigr) - \sum_{x \in A}p(x)\,g\Bigl(\frac{q(x)}{p(x)}\Bigr) \leq \sum_{x \in A}q(x) - 1 \leq 0

である（第 1 の和は絶対収束し、第 2 の和の項は $0$ 以上）。等号なら第 2 の和は $0$ で $\sum_{x \in A}q(x) = 1$ だから、 $A$ 上で $q = p$ 、 $A$ の外で $q = 0 = p$ である。密度の場合も、 $A$ 上で $q = 0$ となる部分は（確率が正なら $\operatorname{KL} = \infty$ なので）零集合として除けば同じ式が成り立ち、等号なら $0$ 以上の関数 $p \cdot g(q/p)$ の $A$ 上の積分が $0$ なので $A$ 上ほとんど至るところ $q = p$ で（06-measure-integration 第3章系 3.12 の 4）、 $\int_A q = 1$ から $A$ の外でもほとんど至るところ $q = 0$ である。 $\square$

系 7.3

$H(p) < \infty$ なら $H(p, q) \geq H(p)$ で、等号は $q = p$ に限る。すなわち $q \mapsto H(p, q)$ は $q = p$ でだけ最小になる。

$K$ 個の元からなる集合の上では $0 \leq H(p) \leq \log K$ で、右の等号は一様分布に限る。

証明. 1 は $H(p, q) = H(p) + \operatorname{KL}(p \Vert q)$ と定理 7.2 による。2 の左は各項が $0$ 以上であることにより、右は一様分布 $u$ について $\operatorname{KL}(p \Vert u) = \log K - H(p)$ となることによる。 $\square$

例 7.4（KL ダイバージェンスの計算）表の確率が $0.5$ , $0.9$ のベルヌーイ分布を $p, q$ とすると $\operatorname{KL}(p \Vert q) = 0.5\log\frac{0.5}{0.9} + 0.5\log\frac{0.5}{0.1} = \log\frac{5}{3} \approx 0.5108$ 、 $\operatorname{KL}(q \Vert p) = 0.9\log 1.8 + 0.1\log 0.2 \approx 0.3681$ で、向きによって値が違う。正規分布どうしでは、密度の比の対数 $\log\frac{\sigma_2}{\sigma_1} - \frac{(x - \mu_1)^2}{2\sigma_1^2} + \frac{(x - \mu_2)^2}{2\sigma_2^2}$ の期待値をとり、 $E[(X - \mu_2)^2] = \sigma_1^2 + (\mu_1 - \mu_2)^2$ を使って

\operatorname{KL}\bigl(N(\mu_1, \sigma_1^2) \Vert N(\mu_2, \sigma_2^2)\bigr) = \log\frac{\sigma_2}{\sigma_1} + \frac{\sigma_1^2 + (\mu_1 - \mu_2)^2}{2\sigma_2^2} - \frac{1}{2} \tag{1}

を得る。例えば $\operatorname{KL}(N(0, 1) \Vert N(1, 4)) = \log 2 - \frac{1}{4} \approx 0.4431$ 、逆向きは $2 - \log 2 \approx 1.3069$ である（数値積分でも確かめた）。共分散行列が対角の多変量正規分布どうしなら、成分ごとの (1) の和になる。KL ダイバージェンスは対称でなく、三角不等式も満たさない（問題 7.1）ので、距離ではない。

定義 7.5（相互情報量）離散確率変数 $X, Y$ の同時分布を $p(x, y)$ とし、 $X, Y$ を独立にした同時分布 $p(x)p(y)$ を $p_X \otimes p_Y$ と書く。 $I(X; Y) = \operatorname{KL}(p \Vert p_X \otimes p_Y)$ を相互情報量 (mutual information)、 $H(X \mid Y) = -\sum_{x, y}p(x, y)\log p(x \mid y)$ を条件付きエントロピーという（和は $p(x, y) > 0$ の組についてとる）。

命題 7.6

$I(X; Y) = I(Y; X) \geq 0$ で、等号は $X$ と $Y$ が独立のときに限る。

$H(X) < \infty$ なら $I(X; Y) = H(X) - H(X \mid Y)$ で、特に $H(X \mid Y) \leq H(X)$ 。

証明. 1 は定理 7.2 と、独立性がすべての $x, y$ で $p(x, y) = p(x)p(y)$ となることと同値であること（22-statistics 第1章命題 1.4）による。2 は $\log\frac{p(x, y)}{p(x)p(y)} = \log p(x \mid y) - \log p(x)$ の期待値をとればよい（ $H(X)$ が有限なので和を分けてよい）。 $\square$

相互情報量は $Y$ を知ることで減る $X$ の不確かさである（問題 7.2）。相関係数と違って、 $0$ になるのは独立なときに限るので、非線形な依存も捉える。例えば $X$ が $\lbrace -1, 0, 1 \rbrace$ 上の一様分布で $Y = X^2$ なら、相関係数は $0$ だが、 $Y$ は $X$ で決まるので $H(Y \mid X) = 0$ で、命題 7.6 を $X, Y$ を入れ替えて使うと $I(X; Y) = H(Y) = \log 3 - \frac{2}{3}\log 2 \approx 0.637 > 0$ である。

7.2 最尤推定と KL ダイバージェンス

定理 7.7（最尤推定と KL ダイバージェンス）データ $x_1, \dots, x_n$ をパラメータ $\theta$ のモデル $p_\theta$ （22-statistics 第3章の $f(x; \theta)$ ）で表し、 $\ell(\theta) = \sum_i\log p_\theta(x_i)$ とする。

データが有限集合または可算集合に値をとるとき、経験分布 $\hat{p}_n(x) = \frac{1}{n}\lvert \lbrace i \mid x_i = x \rbrace \rvert$ について $\frac{1}{n}\ell(\theta) = -H(\hat{p}_n, p_\theta) = -H(\hat{p}_n) - \operatorname{KL}(\hat{p}_n \Vert p_\theta)$ である。したがって、 $\theta$ が対数尤度を最大にすることと $\operatorname{KL}(\hat{p}_n \Vert p_\theta)$ を最小にすることは同値である。

$X \sim p^{\ast}$ で $H(p^{\ast})$ （密度ならその微分エントロピー）が有限なら $E[\log p_\theta(X)] = -H(p^{\ast}) - \operatorname{KL}(p^{\ast} \Vert p_\theta)$ である。したがって対数尤度の期待値の最大化は $\operatorname{KL}(p^{\ast} \Vert p_\theta)$ の最小化と同値で、 $p^{\ast} = p_{\theta_0}$ なら $\theta_0$ は最大点である。

証明. 1：同じ値の項をまとめると $\frac{1}{n}\sum_i\log p_\theta(x_i) = \sum_x\hat{p}_n(x)\log p_\theta(x) = -H(\hat{p}_n, p_\theta)$ 。 $\hat{p}_n$ は高々 $n$ 点に確率をもつので $H(\hat{p}_n) \leq \log n$ （系 7.3）は有限で $\theta$ によらず、7.1 節の分解が使える。2 も同様で、最後の主張は定理 7.2 による。 $\square$

$x_1, \dots, x_n$ が $p^{\ast}$ からの i.i.d. 標本の実現値で $E[\lvert \log p_\theta(X) \rvert] < \infty$ なら、大数の法則（22-statistics 第1章定理 1.27）より $\frac{1}{n}\ell(\theta)$ は $E[\log p_\theta(X)]$ に確率収束するので、最尤推定は真の分布に KL ダイバージェンスの意味で最も近いモデルを標本から探している（22-statistics 第3章定理 3.32 の一致性の証明の概略に現れる不等式は $\operatorname{KL}(p_{\theta_0} \Vert p_\theta) > 0$ そのものである）。真の分布がモデルに含まれなくても、パラメータ空間がコンパクトで、各 $x$ で $\theta \mapsto \log p_\theta(x)$ が連続、 $E[\sup_\theta \lvert \log p_\theta(X) \rvert] < \infty$ であり、 $E[\log p_\theta(X)]$ の最大点（ $H(p^{\ast})$ が有限なら $\operatorname{KL}(p^{\ast} \Vert p_\theta)$ の最小点） $\theta^{\ast}$ がただ一つなら、最尤推定量は $\theta^{\ast}$ に確率収束する（White, 1982 年。主張のみ。一様な大数の法則を使って示される）。連続分布でも、 $-\frac{1}{n}\ell(\theta)$ は経験分布で期待値をとった交差エントロピーである。第1章の言葉では、密度推定は対数損失 $-\log q(x)$ の経験リスク最小化で、そのリスク $H(p^{\ast}, q)$ は（ $H(p^{\ast})$ が有限なら） $q = p^{\ast}$ でだけ最小になる（系 7.3）。分類でも同様である。

命題 7.8（交差エントロピー損失） $\mathcal{Y}$ を有限集合、 $p^{\ast}(\cdot \mid x)$ を $X = x$ のもとでの $Y$ の真の条件付き分布とする。各 $x$ に $\mathcal{Y}$ 上の分布 $q(\cdot \mid x)$ を対応させる予測の、損失 $-\log q(y \mid x)$ についてのリスクは $E[H(p^{\ast}(\cdot \mid X))] + E[\operatorname{KL}(p^{\ast}(\cdot \mid X) \Vert q(\cdot \mid X))]$ であり、 $q(\cdot \mid x) = p^{\ast}(\cdot \mid x)$ となる予測はベイズ最適である。

証明. $X = x$ のもとでの条件付きリスクは $H(p^{\ast}(\cdot \mid x), q(\cdot \mid x))$ で、 $H(p^{\ast}(\cdot \mid x)) \leq \log\lvert \mathcal{Y} \rvert$ は有限だから 7.1 節の分解が使える。第1章の命題 1.6・補題 1.8 と系 7.3 から従う。 $\square$

第2章の命題 2.15 は $\lvert \mathcal{Y} \rvert = 2$ の場合である。0-1 損失（第1章の定理 1.10）と違い、交差エントロピー損失は確率そのものを当てることを求める。

例 7.9（KL ダイバージェンスの向き）最尤推定が最小にする $\operatorname{KL}(p^{\ast} \Vert p_\theta)$ は $p^{\ast}$ で期待値をとるので、データのあるところで $p_\theta \approx 0$ だと大きな罰を受ける。 $p^{\ast} = \frac{1}{2}N(-3, 1) + \frac{1}{2}N(3, 1)$ を $N(m, s^2)$ で近似すると、 $\operatorname{KL}(p^{\ast} \Vert N(m, s^2))$ の最小点は平均と分散を合わせた $N(0, 10)$ で（問題 7.3）、データのほとんどない $0$ の付近に密度の山をおく。逆向きの $\operatorname{KL}(N(m, s^2) \Vert p^{\ast})$ は近似する側で期待値をとるので $p^{\ast} \approx 0$ のところを避け、数値的に最小化すると、最小点は $N(\pm 2.98, 1.05)$ 付近（片方の山だけ。値は約 $0.689$ ）になる。逆向きの KL ダイバージェンスは 7.4 節の変分推論に現れる。

7.3 混合ガウスモデルと EM アルゴリズム

定義 7.10（混合ガウスモデル, Gaussian mixture model）混合係数 $\pi_k > 0$ （ $\sum_{k=1}^{K}\pi_k = 1$ ）、平均 $\mu_k \in \mathbb{R}^d$ 、正定値な共分散行列 $\Sigma_k$ の組を $\theta$ とし、 $N_d(\mu, \Sigma)$ の密度（22-statistics 第1章定理 1.22）を $\varphi(x; \mu, \Sigma)$ として、 $p_\theta(x) = \sum_{k=1}^{K}\pi_k\varphi(x; \mu_k, \Sigma_k)$ を密度とする分布を混合ガウス分布という。

確率 $\pi_k$ で成分の番号 $Z = k$ を選んでから $X \sim N_d(\mu_k, \Sigma_k)$ を生成すると、 $X$ の分布は $p_\theta$ 、 $(X, Z)$ の同時分布は $p_\theta(x, k) = \pi_k\varphi(x; \mu_k, \Sigma_k)$ である。潜在変数 $Z$ の事後確率 $\gamma_k(x) = \pi_k\varphi(x; \mu_k, \Sigma_k)/p_\theta(x)$ を負担率 (responsibility) という。対数尤度 $\ell(\theta) = \sum_i\log\sum_k\pi_k\varphi(x_i; \mu_k, \Sigma_k)$ は閉じた形では最大化できないが、各 $x_i$ の成分の番号がわかっていれば、成分ごとの割合・平均・共分散行列を計算するだけでよい。EM アルゴリズムは、わからない番号を、現在のパラメータでの事後分布で平均して埋めることを繰り返す。

一般に、観測 $x$ と有限集合に値をとる潜在変数 $z$ の同時分布 $p_\theta(x, z) > 0$ をもつモデルを考える（ $x, z$ はデータ全体・潜在変数全体でよい）。 $\ell(\theta) = \log p_\theta(x)$ 、 $p_\theta(x) = \sum_zp_\theta(x, z)$ である。

定義 7.11（EM アルゴリズム, expectation–maximization algorithm）初期値 $\theta^{(0)}$ から、 $t = 0, 1, 2, \dots$ について次を繰り返す。

E ステップ：事後分布 $p_{\theta^{(t)}}(z \mid x)$ を求め、 $Q(\theta \mid \theta^{(t)}) = \sum_zp_{\theta^{(t)}}(z \mid x)\log p_\theta(x, z)$ とおく。

M ステップ： $Q(\theta \mid \theta^{(t)})$ を最大にする $\theta$ を $\theta^{(t+1)}$ とする。

データが i.i.d. なら事後分布も $\prod_ip_\theta(z_i \mid x_i)$ と分かれ、 $Q = \sum_i\sum_{z_i}p_{\theta^{(t)}}(z_i \mid x_i)\log p_\theta(x_i, z_i)$ となる。

定理 7.12（EM アルゴリズムの単調性） $Q(\theta \mid \theta') \geq Q(\theta' \mid \theta')$ ならば $\ell(\theta) \geq \ell(\theta')$ である。特に EM アルゴリズムの列について $\ell(\theta^{(0)}) \leq \ell(\theta^{(1)}) \leq \ell(\theta^{(2)}) \leq \cdots$ 。

証明. $w(z) = p_{\theta'}(z \mid x)$ とおくと $p_{\theta'}(x, z) = p_{\theta'}(x)w(z)$ なので、 $\frac{p_\theta(x)}{p_{\theta'}(x)} = \sum_zw(z)\frac{p_\theta(x, z)}{p_{\theta'}(x, z)}$ である。 $w$ は確率分布で $\log$ は凹だから、イェンセンの不等式（01-calculus 第4章定理 4.31 を凸関数 $-\log$ に使う）より

\ell(\theta) - \ell(\theta') \geq \sum_zw(z)\log\frac{p_\theta(x, z)}{p_{\theta'}(x, z)} = Q(\theta \mid \theta') - Q(\theta' \mid \theta') \geq 0

M ステップの $\theta^{(t+1)}$ は $Q(\theta^{(t+1)} \mid \theta^{(t)}) \geq Q(\theta^{(t)} \mid \theta^{(t)})$ を満たす。 $\square$

定理 7.13（混合ガウスモデルの EM）データ $x_1, \dots, x_n \in \mathbb{R}^d$ のすべてを含む超平面はない（ $d = 1$ なら、すべては等しくない）とする。E ステップは現在のパラメータでの負担率 $\gamma_{ik} = \gamma_k(x_i)$ の計算であり、M ステップの最大点はただ一つで、 $N_k = \sum_i\gamma_{ik}$ として
$\pi_k^{\mathrm{new}} = \frac{N_k}{n}, \qquad \mu_k^{\mathrm{new}} = \frac{1}{N_k}\sum_{i=1}^{n}\gamma_{ik}x_i, \qquad \Sigma_k^{\mathrm{new}} = \frac{1}{N_k}\sum_{i=1}^{n}\gamma_{ik}(x_i - \mu_k^{\mathrm{new}})(x_i - \mu_k^{\mathrm{new}})^{\top}$
である。 $\pi_k^{\mathrm{new}} > 0$ で、 $\Sigma_k^{\mathrm{new}}$ は正定値である。

証明. $P(Z_i = k \mid x_i) = \gamma_{ik} > 0$ で、 $Q = \sum_kN_k\log\pi_k + \sum_k\sum_i\gamma_{ik}\log\varphi(x_i; \mu_k, \Sigma_k)$ は $\pi$ の部分と各 $(\mu_k, \Sigma_k)$ の部分に分かれる。 $\sum_k\gamma_{ik} = 1$ より $\nu_k = N_k/n$ は確率分布で、第 1 項は $-nH(\nu, \pi)$ だから、系 7.3 より $\pi = \nu$ でだけ最大になる。 $\bar{x}_k = \mu_k^{\mathrm{new}}$ 、 $S_k = \Sigma_k^{\mathrm{new}}$ とおく。 $\sum_i\gamma_{ik}(x_i - \bar{x}_k) = 0$ で交差項が消えることと $v^{\top}\Sigma_k^{-1}v = \operatorname{tr}(\Sigma_k^{-1}vv^{\top})$ から、第 2 項の $k$ 番目は定数を除いて

-\frac{N_k}{2}\Bigl(\log\det\Sigma_k + \operatorname{tr}(\Sigma_k^{-1}S_k) + (\bar{x}_k - \mu_k)^{\top}\Sigma_k^{-1}(\bar{x}_k - \mu_k)\Bigr)

で、どの $\Sigma_k$ についても $\mu_k = \bar{x}_k$ でだけ最大になる。仮定より $v \neq 0$ なら $v^{\top}x_i$ はすべては等しくないので $v^{\top}S_kv = \frac{1}{N_k}\sum_i\gamma_{ik}(v^{\top}(x_i - \bar{x}_k))^2 > 0$ 、すなわち $S_k$ は正定値である。正の平方根 $S_k^{1/2}$ （02-linear-algebra 第8章命題 8.15）について $S_k^{1/2}\Sigma_k^{-1}S_k^{1/2}$ の固有値を $\lambda_j > 0$ とすると、 $\operatorname{tr}(\Sigma_k^{-1}S_k) = \sum_j\lambda_j$ 、 $\log\det\Sigma_k = \log\det S_k - \sum_j\log\lambda_j$ だから、 $\lambda - \log\lambda \geq 1$ （定理 7.2 の証明）より $\log\det\Sigma_k + \operatorname{tr}(\Sigma_k^{-1}S_k) \geq \log\det S_k + d$ で、等号はすべての $\lambda_j = 1$ 、すなわち $\Sigma_k = S_k$ のときに限る。 $\square$

例 7.14（1 次元の EM）データ $1, 2, 3, 4, 8, 9, 10$ に $K = 2$ の混合ガウス分布を当てはめる。初期値 $\pi = (0.5, 0.5)$ 、 $\mu = (2, 4)$ 、 $\sigma^2 = (1, 1)$ での負担率は $\gamma_{i1} = 1/(1 + e^{2x_i - 6}) = \sigma(6 - 2x_i)$ （第2章のロジスティック関数）で、 $x_i = 1, 2, 3, 4$ に対して $0.982, 0.881, 0.5, 0.119$ 、 $8, 9, 10$ に対しては $10^{-4}$ 未満である。更新を $t$ 回行った後の値は次のとおりである（Python で計算した）。

$t$	$\ell$	$\pi_1$	$\mu_1$	$\mu_2$	$\sigma_1^2$	$\sigma_2^2$
0	$-49.8194$	$0.5000$	$2.0000$	$4.0000$	$1.0000$	$1.0000$
1	$-17.0402$	$0.3546$	$1.9020$	$7.1447$	$0.7804$	$7.4061$
2	$-16.8994$	$0.3879$	$2.0422$	$7.3412$	$0.8448$	$7.1101$
4	$-16.4717$	$0.4658$	$2.2579$	$7.9255$	$1.0707$	$5.3448$
6	$-15.1382$	$0.5514$	$2.4513$	$8.7697$	$1.2200$	$1.7811$
8	$-14.5510$	$0.5714$	$2.5000$	$9.0000$	$1.2501$	$0.6667$

対数尤度は単調に増え、 $\lbrace 1, 2, 3, 4 \rbrace$ と $\lbrace 8, 9, 10 \rbrace$ に分けたときの割合 $4/7$ 、平均 $2.5, 9$ 、分散 $1.25, 2/3$ にほぼ等しい値に落ち着く。

EM が保証するのは、対数尤度が減らないことだけである。 $\nabla_\theta Q(\theta \mid \theta')$ を $\theta = \theta'$ で評価すると $\sum_zp_{\theta'}(z \mid x)\nabla\log p_{\theta'}(x, z) = \nabla\ell(\theta')$ なので、EM が止まる点は（ $Q$ の最大点が内部にあれば） $\ell$ の停留点である。さらに、上位集合 $\lbrace \theta \mid \ell(\theta) \geq \ell(\theta^{(0)}) \rbrace$ がコンパクトで、 $\ell$ が連続かつ内部で微分可能、列がパラメータ空間の内部にとどまり、 $Q(\theta \mid \theta')$ が $(\theta, \theta')$ について連続ならば、EM の列の集積点はすべて $\ell$ の停留点で、 $\ell(\theta^{(t)})$ はある停留点での値に単調に収束する（Wu, 1983 年。主張のみ。混合ガウスモデルでは、次の命題 7.15 のためにこの上位集合がコンパクトでなく、そのままでは使えない）。しかし停留点は最大点とは限らない。例 7.14 で $\mu_1 = \mu_2$ 、 $\sigma_1^2 = \sigma_2^2$ から始めると、すべての負担率が $\pi_k$ に等しくなり、両成分とも 1 つの正規分布の最尤推定値 $\mu_k = 37/7$ 、 $\sigma_k^2 = 556/49$ に移って以後動かない（ $\ell \approx -18.4339 < -14.5510$ ）。そもそも尤度には最大値がない。

命題 7.15（尤度の非有界性） $d = 1$ 、 $K = 2$ とし、データ $x_1, \dots, x_n$ を任意にとる。 $\pi_1 = \pi_2 = 1/2$ 、 $\mu_1 = x_1$ 、 $\mu_2 = 0$ 、 $\sigma_2^2 = 1$ と固定して $\sigma_1 \to +0$ とすると $\ell(\theta) \to \infty$ 。したがって混合ガウスモデルの尤度は上に有界でなく、最尤推定量は存在しない。

証明. $p_\theta(x_1) \geq \frac{1}{2}\varphi(x_1; x_1, \sigma_1^2) = \frac{1}{2\sqrt{2\pi}\sigma_1}$ で、 $p_\theta(x_i) \geq \frac{1}{2}\varphi(x_i; 0, 1)$ は $\sigma_1$ によらない正の数である。よって $\ell(\theta) \geq -\log(2\sqrt{2\pi}\sigma_1) + \sum_{i \geq 2}\log\frac{\varphi(x_i; 0, 1)}{2} \to \infty$ 。 $\square$

成分が 1 点に縮退すると、その点の密度がいくらでも大きくなるのである（1 つの正規分布では、観測値がすべて等しい場合を除いてこうならない。22-statistics 第3章例 3.13）。例 7.14 のデータで初期値を $\pi = (0.5, 0.5)$ 、 $\mu = (1, 6)$ 、 $\sigma^2 = (0.01, 4)$ とすると、1 回の更新で $\sigma_1^2 \approx 2.9 \times 10^{-20}$ 、 $\ell \approx 3.39$ となり、次の更新で $\sigma_1^2$ が倍精度の計算で $0$ になって破綻する。成分の番号を入れ替えても $p_\theta$ は変わらないので、識別可能性も崩れている（22-statistics 第3章注意 3.34）。なお $\Sigma_k = \varepsilon I$ 、 $\pi_k = 1/K$ と固定すると $\gamma_{ik}$ は $e^{-\lVert x_i - \mu_k \rVert^2/(2\varepsilon)}$ に比例するので、 $\varepsilon \to +0$ で最も近い中心の成分で $1$ 、ほかで $0$ に近づき（最も近い中心がただ一つのとき）、EM は第3章の $k$ 平均法（定義 3.18）に近づく。EM は $k$ 平均法の割り当てを確率に「柔らかく」したものである。

ヒント

実務では EM が収束しても、最尤推定値が得られたとは限らない。(1) 初期値を変えて何回も実行し（ $k$ 平均法の結果を初期値にすることも多い）、対数尤度が最大のものを採る。(2) 分散に下限を設けるなどして縮退を防ぐ。同じ値が繰り返し現れるデータ（定価の商品の購入額、丸めた測定値）では特に起こりやすく、異常に大きい対数尤度はまず縮退を疑う。(3) 成分の数 $K$ は訓練データの対数尤度では選べない（成分を複製すれば $K + 1$ 成分で同じ分布を表せるので、達成できる値は $K$ について下がらない）。検証データでの対数尤度や BIC（22-statistics 第6章定義 6.16。最大対数尤度は存在しないので、分散の下限などの制約のもとで EM が見つけた値で代用する）で選ぶ。

7.4 変分推論と ELBO

ベイズ統計（22-statistics 第7章）では未知の量をまとめて $z$ とし、事後分布 $p(z \mid x) = p(x, z)/p(x)$ を求めるが、周辺尤度 $p(x) = \int p(x, z)\ dz$ はたいてい計算できない。変分推論 (variational inference) は、扱いやすい分布の族 $\mathcal{Q}$ から $\operatorname{KL}(q \Vert p(\cdot \mid x))$ が最小の $q$ を探す。この量は未知の $p(x)$ を含むが、次の恒等式で避けられる。

定理 7.16（ELBO 分解） $p(x) > 0$ とし、 $z$ の分布 $q$ は「 $q(z) > 0$ なら $p(x, z) > 0$ 」と $E_q[\lvert \log q(Z) \rvert] < \infty$ を満たすとする。 $\operatorname{ELBO}(q) = E_q[\log p(x, Z)] - E_q[\log q(Z)]$ とおくと
$\log p(x) = \operatorname{ELBO}(q) + \operatorname{KL}\bigl(q \Vert p(\cdot \mid x)\bigr)$
である（KL ダイバージェンスが $\infty$ なら $\operatorname{ELBO}(q) = -\infty$ と読む）。特に $\operatorname{ELBO}(q) \leq \log p(x)$ で、等号は $q = p(\cdot \mid x)$ のときに限る。

証明. $q(z) > 0$ なら $p(z \mid x) = p(x, z)/p(x) > 0$ で、 $\log p(x) = (\log p(x, z) - \log q(z)) + (\log q(z) - \log p(z \mid x))$ である。 $Z \sim q$ について期待値をとると、右辺の第 2 項は $\operatorname{KL}(q \Vert p(\cdot \mid x)) \in [0, \infty]$ 、第 1 項は $\operatorname{ELBO}(q)$ になる。後半は定理 7.2 による。 $\square$

$p(x)$ を証拠 (evidence) ともいうので、 $\operatorname{ELBO}(q)$ を証拠下界 (evidence lower bound) という。ELBO の最大化は KL ダイバージェンスの最小化と同値で、計算には $p(x, z)$ しか要らない。 $z$ の事前分布 $p_Z$ を使えば $\operatorname{ELBO}(q) = E_q[\log p(x \mid Z)] - \operatorname{KL}(q \Vert p_Z)$ （データへの当てはまりと、事前分布から離れることへの罰）とも書ける。最尤推定とは逆向きの KL ダイバージェンスなので、 $q$ は例 7.9 のように事後分布の確率が小さいところを避ける。

定理 7.17（平均場近似の座標上昇） $z = (z_1, \dots, z_m)$ とし、 $q(z) = \prod_jq_j(z_j)$ の形の分布だけを考える（平均場近似, mean-field approximation）。 $j$ 以外の $q_i$ を固定したとき、 $\operatorname{ELBO}(q)$ を最大にする $q_j$ はただ一つで
$q_j^{\ast}(z_j) = \frac{1}{C_j}\exp\Bigl(E_{-j}\bigl[\log p(x, z_j, Z_{-j})\bigr]\Bigr)$
である。ここで $E_{-j}$ は $i \neq j$ の $Z_i \sim q_i$ （互いに独立）についての期待値、 $C_j$ は正規化定数で、これらはすべて有限とする。

証明. $q$ は積の形なので $E_q[\log q(Z)] = \sum_iE_{q_i}[\log q_i(Z_i)]$ であり、 $E_q[\log p(x, Z)]$ は先に $Z_{-j}$ について期待値をとれば $E_{q_j}[\log q_j^{\ast}(Z_j)] + \log C_j$ に等しい。よって $q_j$ によらない定数 $c_j$ について $\operatorname{ELBO}(q) = E_{q_j}[\log q_j^{\ast}(Z_j) - \log q_j(Z_j)] + c_j = -\operatorname{KL}(q_j \Vert q_j^{\ast}) + c_j$ で、定理 7.2 から従う。 $\square$

$j = 1, \dots, m$ の順にこの更新を繰り返すと（座標上昇変分推論, CAVI）ELBO は減らず、 $\log p(x)$ 以下なので値は収束する。ただし ELBO は一般に凹でなく、結果は初期値に依存しうる。

例 7.18（相関のある正規分布）事後分布が $N_2(\mu, \Sigma)$ 、 $\Lambda = \Sigma^{-1}$ のとき、 $\log p(x, z) = -\frac{1}{2}(z - \mu)^{\top}\Lambda(z - \mu) + c$ のうち $z_1$ を含む項は $z_2$ について 1 次なので、定理 7.17 の期待値は $z_2$ を $m_2 = E_{q_2}[Z_2]$ で置き換えて平方完成すれば求まり、 $q_1^{\ast} = N(m_1, 1/\Lambda_{11})$ 、 $m_1 = \mu_1 - \Lambda_{12}(m_2 - \mu_2)/\Lambda_{11}$ となる（ $q_2$ も同様）。相関係数を $\rho$ とすると、平均は $\mu$ に収束するが（誤差は 1 巡ごとに $\rho^2$ 倍になる）、分散 $1/\Lambda_{11} = (1 - \rho^2)\Sigma_{11}$ は本当の周辺分布の分散 $\Sigma_{11}$ より小さく、 $\rho = 0.9$ なら $0.19$ 倍である（計算機でも確かめた）。

注意

変分推論で得た $q$ から作った信用区間は、例 7.18 のように狭すぎることがある。平均場近似は変数の間の相関を表せず、逆向きの KL ダイバージェンスは事後分布の確率が小さいところを避けるからである。不確かさの評価が重要なら、マルコフ連鎖モンテカルロ法の結果と比べて確かめる。

命題 7.19（EM は ELBO の座標上昇である） $\operatorname{ELBO}(q, \theta) = E_q[\log p_\theta(x, Z)] - E_q[\log q(Z)]$ とおく。

$\theta$ を固定すると、 $q$ についての最大点は $q = p_\theta(\cdot \mid x)$ で、最大値は $\ell(\theta)$ である。

$q = p_{\theta^{(t)}}(\cdot \mid x)$ を固定すると、 $\operatorname{ELBO}(q, \theta) = Q(\theta \mid \theta^{(t)}) - E_q[\log q(Z)]$ である。

したがって E ステップは $q$ について、M ステップは $\theta$ についての ELBO の最大化であり、 $\ell(\theta^{(t+1)}) \geq \operatorname{ELBO}(q, \theta^{(t+1)}) \geq \operatorname{ELBO}(q, \theta^{(t)}) = \ell(\theta^{(t)})$ となる。

証明. 1 は定理 7.16 を $p_\theta$ に使えばよく、2 は定義そのものである。最後の不等式は順に 1、M ステップ、1 による。 $\square$

つまり EM は、潜在変数の分布 $q$ に制約を置かず（E ステップで事後分布を厳密に求め）、 $\theta$ は点推定する変分推論である。事後分布が計算できなければ $q$ を平均場近似などに制限し（変分 EM）、 $\theta$ にも事前分布を置いて $z$ に含めれば、ベイズ推論としての変分推論になる。

7.5 隠れマルコフモデル

設備の状態（正常・異常）は直接は見えず、アラームの有無だけが記録される。混合モデルの成分の番号を、時間とともにマルコフ連鎖で変化させたのが隠れマルコフモデルである（音声や品詞の推定などにも使われる）。

定義 7.20（隠れマルコフモデル, hidden Markov model）状態 $\lbrace 1, \dots, K \rbrace$ の初期分布 $\pi$ 、推移確率 $a_{jk} \geq 0$ （ $\sum_ka_{jk} = 1$ ）、各状態 $k$ の出力分布 $b_k$ （確率関数または密度）を与え、状態の列 $z = (z_1, \dots, z_T)$ と観測の列 $x = (x_1, \dots, x_T)$ の同時分布を
$p(x, z) = \pi_{z_1}b_{z_1}(x_1)\prod_{t=2}^{T}a_{z_{t-1}z_t}b_{z_t}(x_t)$
で定めるモデルを隠れマルコフモデル（HMM）という。

状態の列はマルコフ連鎖（11-probability 第6章定義 6.1）で、状態が与えられると観測は独立であり、 $x_t$ は $z_t$ だけに依存する。 $x_{1:t} = (x_1, \dots, x_t)$ と書く。最後の因子 $a_{z_{T-1}z_T}b_{z_T}(x_T)$ を $x_T, z_T$ について和（積分）をとると $1$ になることを繰り返せば、 $t \leq T$ について

p(x_{1:t}, z_{1:t}) = \pi_{z_1}b_{z_1}(x_1)\prod_{s=2}^{t}a_{z_{s-1}z_s}b_{z_s}(x_s) \tag{2}

である。尤度 $p(x) = \sum_zp(x, z)$ は $K^T$ 項の和で、 $K = 10$ , $T = 100$ なら $10^{100}$ 項になる。

定理 7.21（前向きアルゴリズム, forward algorithm） $\alpha_t(k) = p(x_{1:t}, Z_t = k)$ （(2) を $z_{1:t-1}$ について和をとり $z_t = k$ としたもの）とおくと
$\alpha_1(k) = \pi_kb_k(x_1), \qquad \alpha_{t+1}(k) = b_k(x_{t+1})\sum_{j=1}^{K}\alpha_t(j)a_{jk} \quad (1 \leq t < T), \qquad p(x) = \sum_{k=1}^{K}\alpha_T(k)$
であり、 $p(x)$ は $O(TK^2)$ 回の四則演算で計算できる。

証明. (2) より $p(x_{1:t+1}, z_{1:t+1}) = p(x_{1:t}, z_{1:t})a_{z_tz_{t+1}}b_{z_{t+1}}(x_{t+1})$ で、 $z_{t+1} = k$ を固定して $z_{1:t-1}$ と $z_t = j$ について和をとれば漸化式を得る。 $\alpha_1$ は (2) の $t = 1$ の場合で、 $\sum_k\alpha_T(k) = \sum_zp(x, z) = p(x)$ である。各 $t$ で $K$ 個の $\alpha_{t+1}(k)$ をそれぞれ $K + 1$ 回の乗算と $K - 1$ 回の加算で求めるので 1 段は $O(K^2)$ 、全体で $O(TK^2)$ である。 $\square$

素朴な和は、 $K^T$ 個の項のそれぞれに $2T - 1$ 回の乗算が要るので $O(TK^T)$ である。マルコフ性により、先の計算に要る過去の情報が $\alpha_t$ の $K$ 個の数に集約されるのであり、動的計画法（23-optimization 第7章 7.8 節）の一例である。

例 7.22（アラームの記録）状態 1 を正常、2 を異常、観測 $1$ をアラームあり、 $0$ をなしとし、 $\pi = (0.9, 0.1)$ 、 $a_{11} = 0.9$ , $a_{12} = 0.1$ , $a_{21} = 0.3$ , $a_{22} = 0.7$ 、 $b_1(1) = 0.1$ , $b_2(1) = 0.6$ とする。観測 $x = (0, 1, 1)$ では次のようになる。

$t$	$x_t$	$\alpha_t(1)$	$\alpha_t(2)$
1	0	$0.9 \times 0.9 = 0.81$	$0.1 \times 0.4 = 0.04$
2	1	$0.1 \times (0.81 \times 0.9 + 0.04 \times 0.3) = 0.0741$	$0.6 \times (0.81 \times 0.1 + 0.04 \times 0.7) = 0.0654$
3	1	$0.1 \times (0.0741 \times 0.9 + 0.0654 \times 0.3) = 0.008631$	$0.6 \times (0.0741 \times 0.1 + 0.0654 \times 0.7) = 0.031914$

$p(x) = 0.040545$ で、 $8$ 通りの状態の列について $p(x, z)$ を足したものと一致する（計算機で確かめた）。正規化した $\alpha_t(k)/\sum_j\alpha_t(j) = P(Z_t = k \mid x_{1:t})$ （フィルタリング）によると、異常の確率は時刻 2 で $0.0654/0.1395 \approx 0.469$ 、時刻 3 で $0.031914/0.040545 \approx 0.787$ である。

長い系列では $\alpha_t$ が指数的に小さくなる。このモデルで $0$ と $1$ を交互に並べた長さ 1000 の観測列では $\log p(x) \approx -919.65$ 、すなわち $p(x) \approx 10^{-399.4}$ で、倍精度で表せる最小の正の数（約 $4.9 \times 10^{-324}$ ）より小さく、漸化式をそのまま計算すると途中（ $t = 810$ ）で $\alpha_t$ が $0$ になる（計算機で確かめた）。そこで各段で $\alpha_t$ を正規化し、正規化定数の対数を足していく。 $t$ 段目の正規化定数は $p(x_t \mid x_{1:t-1})$ であり、7.6 節の自己回帰分解 $\log p(x) = \sum_t\log p(x_t \mid x_{1:t-1})$ を計算していることになる。

ビタビ・アルゴリズム (Viterbi algorithm)：最も確からしい状態の列 $\arg\max_zp(x, z)$ は、漸化式の和を最大値に変えた $\delta_1(k) = \pi_kb_k(x_1)$ 、 $\delta_{t+1}(k) = b_k(x_{t+1})\max_j\delta_t(j)a_{jk}$ を計算し、最大を与えた $j$ を記録して最後から逆にたどれば、 $O(TK^2)$ で求まる。例 7.22 では（正常, 異常, 異常）で、その事後確率は $0.020412/0.040545 \approx 0.503$ である。

バウム–ウェルチ・アルゴリズム (Baum–Welch algorithm)：パラメータが未知なら EM アルゴリズムで推定する。E ステップに要る $P(Z_t = k \mid x) = \alpha_t(k)\beta_t(k)/p(x)$ などは、後ろ向きの量 $\beta_t(k) = p(x_{t+1:T} \mid Z_t = k)$ を同様の漸化式で求めれば計算でき（問題 7.5）、M ステップでは推移の回数の期待値の比などでパラメータを更新する（導出は省略する）。定理 7.12 により尤度は単調に増加するが（ $a_{jk} = 0$ などで $p_\theta(x, z) = 0$ となる $z$ があっても、その証明の和を $p_{\theta'}(z \mid x) > 0$ の $z$ に限れば、最初の等式が $\geq$ になるだけで成り立つ）、局所解の問題は混合ガウスモデルと同じである。

7.6 生成モデル（紹介）

代表的な 3 つの生成モデルを、定式化と学習の目的関数に絞って本章の言葉で読む（ネットワークの構造や性能には立ち入らない）。

変分オートエンコーダ

変分オートエンコーダ (variational autoencoder, VAE) は、 $z \sim N_m(0, I)$ とニューラルネットワーク $f_\theta$ による $p_\theta(x \mid z) = N_d(f_\theta(z), \sigma^2I)$ などのモデル（デコーダ）を、事後分布の近似 $q_\phi(z \mid x) = N_m(\mu_\phi(x), \operatorname{diag}(s_\phi(x)^2))$ （エンコーダ。これもネットワーク）とともに、ELBO

E_{q_\phi(z \mid x)}\bigl[\log p_\theta(x \mid Z)\bigr] - \operatorname{KL}\bigl(q_\phi(\cdot \mid x) \Vert N_m(0, I)\bigr) \leq \log p_\theta(x)

のデータについての和を最大にして学習する（変分 EM と同じ ELBO を、 $q$ をネットワークで表し、E ステップと M ステップを交互に解く代わりに $\theta$ と $\phi$ について同時に勾配法で最大化する）。第 1 項は定数を除いて $-\frac{1}{2\sigma^2}E\lVert x - f_\theta(Z) \rVert^2$ （ $x$ を $z$ に符号化して復元したときの誤差）、KL ダイバージェンスの項は (1) を成分ごとに足して $\frac{1}{2}\sum_j(\mu_j^2 + s_j^2 - 1 - \log s_j^2)$ （ $\mu = \mu_\phi(x)$ 、 $s = s_\phi(x)$ ）である。学習後は $z \sim N_m(0, I)$ から $f_\theta(z)$ を作れば新しいデータが得られる。期待値をとる分布が $\phi$ によるので、勾配は次の書き換えで計算する。

命題 7.23（再パラメータ化, reparameterization） $\varepsilon \sim N_m(0, I)$ 、 $s_j > 0$ なら、成分ごとの積 $\odot$ について $\mu + s \odot \varepsilon \sim N_m(\mu, \operatorname{diag}(s^2))$ であり、 $E_{Z \sim N_m(\mu, \operatorname{diag}(s^2))}[g(Z)] = E[g(\mu + s \odot \varepsilon)]$ である。 $g$ が $C^1$ 級で微分と期待値の順序を交換できるなら
$\frac{\partial}{\partial \mu_j}E[g(\mu + s \odot \varepsilon)] = E\left[\frac{\partial g}{\partial z_j}(\mu + s \odot \varepsilon)\right], \qquad \frac{\partial}{\partial s_j}E[g(\mu + s \odot \varepsilon)] = E\left[\varepsilon_j\frac{\partial g}{\partial z_j}(\mu + s \odot \varepsilon)\right]$

証明. $\mu + s \odot \varepsilon = \mu + \operatorname{diag}(s)\varepsilon$ だから、前半は 22-statistics 第1章定理 1.22 の 2 による。後半は期待値の中を連鎖律で微分すればよい（交換の十分条件は 06-measure-integration 第3章定理 3.25）。 $\square$

右辺は $\varepsilon$ を生成すれば不偏に推定でき、 $\mu = \mu_\phi(x)$ 、 $s = s_\phi(x)$ を通して誤差逆伝播法（第5章定理 5.5）で $\phi$ の勾配が得られる。

拡散モデル

拡散モデル (diffusion model) は、データ $x_0 \sim p^{\ast}$ に雑音を少しずつ加える前向き過程 $x_t = \sqrt{1 - \beta_t}x_{t-1} + \sqrt{\beta_t}\varepsilon_t$ （ $t = 1, \dots, T$ 。 $0 < \beta_t < 1$ は定数、 $\varepsilon_t$ は $x_0$ と独立な $N_d(0, I)$ の i.i.d.）と、それを逆にたどる、学習する逆過程からなる。

命題 7.24（前向き過程の周辺分布） $\bar{\alpha}_t = \prod_{s=1}^{t}(1 - \beta_s)$ とおくと、 $x_0$ を与えたもとで $x_t \sim N_d(\sqrt{\bar{\alpha}_t}x_0, (1 - \bar{\alpha}_t)I)$ である。すなわち $x_t$ は、 $x_0$ と独立な $\varepsilon \sim N_d(0, I)$ による $\sqrt{\bar{\alpha}_t}x_0 + \sqrt{1 - \bar{\alpha}_t}\varepsilon$ と同じ分布に従う。

証明. $u_t = x_t - \sqrt{\bar{\alpha}_t}x_0$ とおくと、 $\sqrt{1 - \beta_t}\sqrt{\bar{\alpha}_{t-1}} = \sqrt{\bar{\alpha}_t}$ より $u_0 = 0$ 、 $u_t = \sqrt{1 - \beta_t}u_{t-1} + \sqrt{\beta_t}\varepsilon_t$ である。帰納法により $u_t$ は $\varepsilon_1, \dots, \varepsilon_t$ の 1 次結合なので、 $x_0$ と独立で、平均 $0$ の多変量正規分布に従う（22-statistics 第1章定義 1.21・定理 1.22）。共分散行列を $v_tI$ とすると、独立性から $v_t = (1 - \beta_t)v_{t-1} + \beta_t$ 、 $v_0 = 0$ で、 $v_t = 1 - \bar{\alpha}_t$ がこれを満たす。 $\square$

$\bar{\alpha}_T \approx 0$ なら $x_T$ はデータによらずほぼ $N_d(0, I)$ に従う（ $\beta_t = 0.02$ 、 $T = 200$ なら $\bar{\alpha}_T \approx 0.018$ ）。逆過程は $x_T \sim N_d(0, I)$ から正規分布 $p_\theta(x_{t-1} \mid x_t) = N_d(m_\theta(x_t, t), \sigma_t^2I)$ で順に $x_0$ まで生成する（以下、分散 $\sigma_t^2$ は固定する）。 $x_{1:T}$ を潜在変数、前向き過程を学習しない変分分布とみると定理 7.16 の下界が得られ、その符号を変えたものは、 $p_\theta(x_{t-1} \mid x_t)$ と、閉じた形で書ける正規分布

q(x_{t-1} \mid x_t, x_0) = N_d\Bigl(\frac{1}{\sqrt{1 - \beta_t}}\Bigl(x_t - \frac{\beta_t}{\sqrt{1 - \bar{\alpha}_t}}\varepsilon\Bigr), \frac{(1 - \bar{\alpha}_{t-1})\beta_t}{1 - \bar{\alpha}_t}I\Bigr) \qquad (2 \leq t \leq T)

との KL ダイバージェンスの和などに分かれる。ここで $\varepsilon = (x_t - \sqrt{\bar{\alpha}_t}x_0)/\sqrt{1 - \bar{\alpha}_t}$ であり、この式は、 $q(x_t \mid x_{t-1})$ と命題 7.24 の $q(x_{t-1} \mid x_0)$ の積を $x_{t-1}$ について平方完成すれば得られる。平均 $m_\theta(x_t, t)$ を、この式の平均に現れる $\varepsilon$ を雑音を予測するネットワーク $\varepsilon_\theta(x_t, t)$ に置き換えたものにとると、(1) を成分ごとに使えば各項は $\frac{\beta_t^2}{2\sigma_t^2(1 - \beta_t)(1 - \bar{\alpha}_t)}E\lVert \varepsilon - \varepsilon_\theta(\sqrt{\bar{\alpha}_t}x_0 + \sqrt{1 - \bar{\alpha}_t}\varepsilon, t) \rVert^2$ に $\theta$ によらない定数を加えたものになる（ $t = 1$ の項 $-E[\log p_\theta(x_0 \mid x_1)]$ も同じ形で、 $\operatorname{KL}(q(x_T \mid x_0) \Vert N_d(0, I))$ は $\theta$ によらない）。実際には、 $t$ を一様に選び、この係数を省いた二乗誤差を最小にすることが多い（係数を省くのは経験的な選択である）。第1章の定理 1.9 より最適な予測は $E[\varepsilon \mid x_t]$ で、これは $x_t$ の密度 $q_t$ の対数の勾配（スコア）を使って $-\sqrt{1 - \bar{\alpha}_t}\nabla\log q_t(x_t)$ と表せる（問題 7.6）。刻み幅 $h$ について $\beta_t = \beta(th)h$ とおいて $h \to 0$ とする極限では、前向き過程は確率微分方程式 $dX_s = -\frac{1}{2}\beta(s)X_s\ ds + \sqrt{\beta(s)}\ dB_s$ の解に近づく（主張のみ）。 $\beta$ が定数なら、各成分は 11-probability 第7章例 7.20 のオルンシュタイン–ウーレンベック過程で、同例より、出発点によらず分散 $(\sqrt{\beta})^2/(2 \cdot \beta/2) = 1$ の正規分布 $N(0, 1)$ に近づく。

大規模言語モデル

文章を、有限集合 $V$ （語彙）の元であるトークン（単語やその断片）の列 $x_1, \dots, x_T$ で表し、 $x_{< t} = (x_1, \dots, x_{t-1})$ と書く。

命題 7.25（自己回帰分解） $V^T$ 上の任意の確率分布 $p$ と、 $p(x_1, \dots, x_T) > 0$ となる列について
$p(x_1, \dots, x_T) = \prod_{t=1}^{T}p(x_t \mid x_{< t})$
が成り立つ。ここで $p(x_t \mid x_{< t}) = p(x_{\leq t})/p(x_{< t})$ は最初の $t$ 個と $t - 1$ 個の周辺分布の比である（ $t = 1$ では $p(x_1)$ ）。

証明. $p(x_{\leq t}) \geq p(x_1, \dots, x_T) > 0$ なので比が定義でき、積をとると隣り合う分母と分子が打ち消し合って $p(x_{\leq T})$ が残る。 $\square$

この分解には何の仮定も要らない。自己回帰型の言語モデルは、各条件付き分布を、前のトークンからニューラルネットワーク（現在の大規模言語モデルの多くはトランスフォーマー。第5章 5.9 節）で計算した特徴量 $h_\theta(x_{< t})$ のソフトマックス回帰 $p_\theta(x_t \mid x_{< t}) = \operatorname{softmax}(Wh_\theta(x_{< t}))_{x_t}$ （第2章の定義 2.19）で表し、大規模言語モデルはこれを非常に多くのパラメータと大量の文章で学習したものである。学習で最小にする平均の交差エントロピー $L(\theta) = -\frac{1}{N}\sum_{t=1}^{N}\log p_\theta(x_t \mid x_{< t})$ （ $N$ はトークン数）は、命題 7.25 より文章全体の負の対数尤度の $1/N$ なので、これは最尤推定である（定理 7.7・命題 7.8）。

定義 7.26（パープレキシティ, perplexity）評価用の文章の $N$ 個のトークンについて、 $\operatorname{PPL} = \exp\bigl(-\frac{1}{N}\sum_{t=1}^{N}\log p_\theta(x_t \mid x_{< t})\bigr)$ をパープレキシティという。

パープレキシティは、正解のトークンに与えた確率の逆数の幾何平均で、 $1$ 以上である。すべてのトークンに確率 $1/\lvert V \rvert$ を与えるモデルではちょうど $\lvert V \rvert$ なので、「平均して何個の候補の間で迷っているか」と読める。文章の生成は $x_t \sim p_\theta(\cdot \mid x_{< t})$ を順に生成して行う（温度 $\tau > 0$ を使って $\operatorname{softmax}(z/\tau)$ から生成すると、 $\tau$ が小さいほど確率の高いトークンに集中する）。対話に使うモデルでは、この最尤推定（事前学習）の後に、対話の例についての同じ交差エントロピーの最小化や、人の評価を使う別の目的関数などで追加の学習をするのが一般的である。

注意

最尤推定が求めるのは学習データの文章の分布をまねることであり、生成された文章の内容が事実として正しいことは、目的関数に直接は入っていない。パープレキシティが低いことも、内容の正しさを意味しない。

ヒント

実務では パープレキシティはトークンあたりの量なので、トークンへの分け方が違うモデルどうしでは比べられない（細かく分けるほど 1 トークンあたりの予測は易しくなる）。比べるなら、同じ評価用の文章全体の対数尤度を文字数やバイト数で割った量を使う（問題 7.7）。評価用の文章が学習データに含まれていると評価は楽観的になる（第1章のデータ漏洩）。長い文章の確率はすぐに浮動小数点数の範囲を下回るので、計算はつねに対数で行う。

まとめ

$H(p, q) = H(p) + \operatorname{KL}(p \Vert q)$ であり、ギブスの不等式 $\operatorname{KL}(p \Vert q) \geq 0$ （等号は $p = q$ ）は $\log u \leq u - 1$ から従う。相互情報量は、独立な場合の分布からの KL ダイバージェンスである。
最尤推定は経験分布との KL ダイバージェンス（交差エントロピー）の最小化であり、交差エントロピー損失のベイズ最適な予測は真の条件付き分布である。
混合ガウスモデルの EM アルゴリズムは負担率の計算（E）と重みつき最尤推定（M）の繰り返しで、対数尤度は単調に増加する。しかし最大でない停留点で止まることがあり、尤度は非有界で、結果は初期値に依存する。
$\log p(x) = \operatorname{ELBO}(q) + \operatorname{KL}(q \Vert p(\cdot \mid x))$ 。平均場近似の更新は $q_j \propto \exp(E_{-j}[\log p(x, Z)])$ で、広がりを過小評価しやすい。EM は $q$ に制約を置かない ELBO の座標上昇である。
隠れマルコフモデルの尤度は前向きアルゴリズムで $O(TK^2)$ で計算でき（素朴な和は $O(TK^T)$ ）、状態の列の推定はビタビ・アルゴリズム、パラメータの推定はバウム–ウェルチ・アルゴリズム（EM）で行う。
変分オートエンコーダは ELBO と再パラメータ化、拡散モデルは閉じた形の前向き過程と雑音の予測、大規模言語モデルは自己回帰分解と交差エントロピー（最尤推定）で学習する。パープレキシティは平均交差エントロピーの指数である。

演習問題

問題 7.1 ★ (1) $p = (1/2, 1/4, 1/4)$ と一様分布 $u = (1/3, 1/3, 1/3)$ について、 $H(p)$ 、 $\operatorname{KL}(p \Vert u)$ 、 $\operatorname{KL}(u \Vert p)$ を求めよ。(2) 表の確率が $0.1, 0.5, 0.9$ のベルヌーイ分布を $p_1, p_2, p_3$ とする。 $\operatorname{KL}(p_1 \Vert p_3) > \operatorname{KL}(p_1 \Vert p_2) + \operatorname{KL}(p_2 \Vert p_3)$ を確かめよ。

解答

(1) $H(p) = \frac{1}{2}\log 2 + 2 \cdot \frac{1}{4}\log 4 = \frac{3}{2}\log 2 \approx 1.0397$ 、 $\operatorname{KL}(p \Vert u) = \log 3 - H(p) \approx 0.0589$ （系 7.3 の証明）、 $\operatorname{KL}(u \Vert p) = \frac{1}{3}(\log\frac{2}{3} + 2\log\frac{4}{3}) = \frac{1}{3}\log\frac{32}{27} \approx 0.0566$ 。

(2) $\operatorname{KL}(p_1 \Vert p_3) = 0.1\log\frac{0.1}{0.9} + 0.9\log\frac{0.9}{0.1} = 0.8\log 9 \approx 1.7578$ 。一方 $\operatorname{KL}(p_1 \Vert p_2) = 0.1\log 0.2 + 0.9\log 1.8 \approx 0.3681$ 、 $\operatorname{KL}(p_2 \Vert p_3) = \log\frac{5}{3} \approx 0.5108$ （例 7.4）で、和は約 $0.8789$ にすぎない。

問題 7.2 ★ メールに特定の単語が含まれるかを $X$ （含めば $1$ ）、迷惑メールかを $Y$ （迷惑メールなら $1$ ）とし、 $P(X = 1, Y = 1) = P(X = 0, Y = 0) = 0.4$ 、 $P(X = 1, Y = 0) = P(X = 0, Y = 1) = 0.1$ とする。 $I(X; Y)$ 、 $H(X)$ 、 $H(X \mid Y)$ を求め、命題 7.6 の 2 を確かめよ。

解答

周辺分布はどちらも $(1/2, 1/2)$ なので $H(X) = \log 2 \approx 0.6931$ 、 $I(X; Y) = 0.8\log\frac{0.4}{0.25} + 0.2\log\frac{0.1}{0.25} = 0.8\log 1.6 + 0.2\log 0.4 \approx 0.1927$ 。 $Y$ のどちらの値のもとでも $X$ の条件付き分布は $(0.8, 0.2)$ なので、 $H(X \mid Y) = -0.8\log 0.8 - 0.2\log 0.2 \approx 0.5004$ で、 $H(X) - H(X \mid Y) \approx 0.1927 = I(X; Y)$ である。

問題 7.3 ★★ $p$ を平均 $m$ 、分散 $v > 0$ の $\mathbb{R}$ 上の密度で、微分エントロピー $h(p) = -\int p\log p$ が有限なものとする。 $N(\mu, s^2)$ のうち $\operatorname{KL}(p \Vert N(\mu, s^2))$ を最小にするのは $\mu = m$ 、 $s^2 = v$ であることを示し、例 7.9 の $N(0, 10)$ を確かめよ。

解答

$-\log\varphi(x; \mu, s^2) = \frac{1}{2}\log(2\pi s^2) + \frac{(x - \mu)^2}{2s^2}$ と $E_p[(X - \mu)^2] = v + (m - \mu)^2$ より

\operatorname{KL}(p \Vert N(\mu, s^2)) = -h(p) + \frac{1}{2}\log(2\pi s^2) + \frac{v + (m - \mu)^2}{2s^2}

で、これは $\mu = m$ でだけ最小になる。そのとき $u = s^2$ の関数 $\frac{1}{2}\log u + \frac{v}{2u}$ の導関数 $\frac{u - v}{2u^2}$ は $u = v$ の前後で負から正に変わる。例 7.9 の混合分布の平均は $0$ 、分散は $E[X^2] = 1 + 9 = 10$ である。

問題 7.4 ★★ ある分析者が、顧客の購入額の対数に混合ガウスモデルを当てはめ、成分の数 $K = 1, \dots, 8$ のそれぞれで EM アルゴリズムを 1 回ずつ実行した。訓練データの対数尤度は $K$ とともにおおむね増え、 $K = 8$ でほかより格段に大きくなった。 $K = 8$ の解では、ある成分の分散が $10^{-9}$ で、その重みは $1/n$ に近かった。分析者は「顧客は 8 つの層に分かれる」と結論した。この分析の問題点を挙げ、どうすべきかを述べよ。

解答

(1) 訓練データの対数尤度は、成分を複製すれば $K$ について下がらないので、 $K$ の選択に使えない。(2) 分散が $10^{-9}$ で重みが $1/n$ に近い成分はほぼ 1 点に縮退しており、格段に大きい対数尤度は命題 7.15 の特異点による見かけのものである（同じ購入額が繰り返し現れると起こりやすい）。(3) 各 $K$ で 1 回しか実行しておらず、局所解や初期値への依存を確かめていない。(4) 成分は顧客の層とは限らない。山が 1 つでも、歪んだ分布や裾の重い分布はいくつもの正規分布の和で近似される。分散に下限を設け、初期値を変えて何回も実行し、 $K$ は検証データでの対数尤度や BIC で選び、解の安定性と、成分が業務上意味のある違いに対応するかを確かめるべきである。

問題 7.5 ★★ 隠れマルコフモデルで $\beta_T(k) = 1$ 、 $t < T$ について $\beta_t(k) = \sum_{z_{t+1}, \dots, z_T}\prod_{s=t+1}^{T}a_{z_{s-1}z_s}b_{z_s}(x_s)$ （ $z_t = k$ とする）と定める。(1) $\beta_t(j) = \sum_ka_{jk}b_k(x_{t+1})\beta_{t+1}(k)$ と $p(x, Z_t = k) = \alpha_t(k)\beta_t(k)$ を示せ。(2) 例 7.22 で $\beta_1, \beta_2$ と $P(Z_2 = 2 \mid x)$ を求め、時刻 2 のフィルタリングの値と比べよ。

解答

(1) $\beta_t(j)$ の定義の和で $z_{t+1} = k$ を固定すると因子 $a_{jk}b_k(x_{t+1})$ が外に出て、残りの和が $\beta_{t+1}(k)$ になる。 $z_t = k$ のとき $p(x, z)$ は $s \leq t$ の因子の積（(2) の右辺）と $s > t$ の因子の積の積であり、 $z_{1:t-1}$ について和をとると前者は $\alpha_t(k)$ 、 $z_{t+1:T}$ について和をとると後者は $\beta_t(k)$ になる。

(2) $\beta_2(1) = 0.9 \times 0.1 + 0.1 \times 0.6 = 0.15$ 、 $\beta_2(2) = 0.3 \times 0.1 + 0.7 \times 0.6 = 0.45$ 、 $\beta_1(1) = 0.9 \times 0.1 \times 0.15 + 0.1 \times 0.6 \times 0.45 = 0.0405$ 、 $\beta_1(2) = 0.3 \times 0.1 \times 0.15 + 0.7 \times 0.6 \times 0.45 = 0.1935$ （検算： $0.81 \times 0.0405 + 0.04 \times 0.1935 = 0.040545 = p(x)$ ）。 $P(Z_2 = 2 \mid x) = 0.0654 \times 0.45/0.040545 \approx 0.726$ で、時刻 2 までの観測によるフィルタリングの値 $0.469$ より大きい。時刻 3 のアラームが、時刻 2 の異常の証拠を強めたのである。

問題 7.6 ★★★ $X_0$ を密度 $p_0$ をもつ $\mathbb{R}^d$ の確率ベクトル、 $\varepsilon \sim N_d(0, I)$ を $X_0$ と独立とし、 $0 < \bar{\alpha} < 1$ について $X = \sqrt{\bar{\alpha}}X_0 + \sqrt{1 - \bar{\alpha}}\varepsilon$ の密度を $q$ とする。微分と積分の順序交換を認めて、 $E[\varepsilon \mid X = x] = -\sqrt{1 - \bar{\alpha}}\nabla\log q(x)$ （ツイーディーの公式）を示せ。

解答

$X_0 = x_0$ のもとでの $X$ の密度 $k(x \mid x_0) = \varphi(x; \sqrt{\bar{\alpha}}x_0, (1 - \bar{\alpha})I)$ について $q(x) = \int k(x \mid x_0)p_0(x_0)\ dx_0$ 、 $\nabla_xk(x \mid x_0) = -\frac{x - \sqrt{\bar{\alpha}}x_0}{1 - \bar{\alpha}}k(x \mid x_0)$ である。 $k(x \mid x_0)p_0(x_0)/q(x)$ は $X = x$ のもとでの $X_0$ の条件付き密度なので

\nabla\log q(x) = \frac{\nabla q(x)}{q(x)} = -\frac{1}{1 - \bar{\alpha}}\int(x - \sqrt{\bar{\alpha}}\,x_0)\frac{k(x \mid x_0)p_0(x_0)}{q(x)}\ dx_0 = -\frac{1}{1 - \bar{\alpha}}E\bigl[X - \sqrt{\bar{\alpha}}\,X_0 \bigm| X = x\bigr]

で、 $X - \sqrt{\bar{\alpha}}X_0 = \sqrt{1 - \bar{\alpha}}\varepsilon$ を代入すればよい。雑音を予測するネットワークは、雑音を加えたデータの分布のスコアを推定しているのである。

問題 7.7 ★★ (1) ある言語モデルが評価用の文章の 4 個のトークンに与えた確率が順に $1/2, 1/4, 1/8, 1/2$ だった。平均交差エントロピーとパープレキシティを求めよ。(2) 同じ評価用の文章で、語彙が 32000 のモデル A のパープレキシティは 20、1 文字を 1 トークンとするモデル B のパープレキシティは 8 だった。A の 1 トークンは平均 1.5 文字にあたる。「B のほうが文章をよく予測している」という結論は正しいか。

解答

(1) 平均交差エントロピーは $\frac{1}{4}(\log 2 + \log 4 + \log 8 + \log 2) = \frac{7}{4}\log 2 \approx 1.213$ 、パープレキシティは $2^{7/4} \approx 3.364$ 。

(2) 正しくない。文章の文字数を $C$ とすると、文章全体の対数尤度は A では $-\frac{C}{1.5}\log 20$ 、B では $-C\log 8$ で、1 文字あたり A は約 $-1.997$ 、B は約 $-2.079$ である。A のほうが文章全体に高い確率を与えており（1 文字あたりのパープレキシティは $20^{2/3} \approx 7.37 < 8$ ）、よく予測している。トークンへの分け方が違うモデルを、トークンあたりのパープレキシティで比べてはいけない。