第7章ベイズ統計

この章の目標

• 事前分布と尤度から事後分布を求め、ベータ–二項・ガンマ–ポアソン・正規–正規の共役事前分布の計算を証明できる
• 二乗損失のベイズ推定量が事後平均、絶対値損失のベイズ推定量が事後中央値であることを証明できる
• 信用区間と信頼区間の解釈の違いを説明し、事後予測分布を計算できる
• メトロポリス–ヘイスティングス法が目標分布を定常分布にもつことを詳細つり合いから証明し、収束にはさらに条件が要ることを説明できる
• ギブスサンプリング、階層モデルによる縮小推定、ベイズファクターの考え方と注意点を説明できる

前提：第1章、第3章、第4章。7.6 節ではマルコフ連鎖の定常分布と収束の定理を 11 確率論 第6章 から引用する。7.8 節では 第6章 の BIC と比べる。

新しい商品ページを公開したところ、20 人が訪れて 3 人が購入した。購入率の最尤推定値は $3/20 = 0.15$ だが（第3章）、20 人では心もとない。一方、同じサイトの過去のページの購入率がおおむね 5〜20% だったことはわかっている。この知識を推定に組み込みたい。また、意思決定をする人が知りたいのは「購入率が 10% を超える確率」のような、パラメータについての確率であることが多い。

これまでの章の方法（頻度論, frequentist）では、パラメータ $\theta$ は未知の定数であり、確率はデータの揺らぎを表すためだけに使った。ベイズ統計 (Bayesian statistics) では、$\theta$ についての不確かさも確率分布で表し、データを見る前の分布（事前分布）を、ベイズの定理によってデータを見た後の分布（事後分布）に更新する。推定・区間・予測は、すべて事後分布から導かれる。本章では、計算が閉じた形でできる共役事前分布、ベイズ推定量、信用区間、事後予測分布を扱った後、一般の事後分布に従う乱数を作るマルコフ連鎖モンテカルロ法が正しい分布を目標にしていることを証明する。

7.1 事前分布と事後分布

本章ではパラメータ $\theta \in \Theta \subset \mathbb{R}^k$ も確率変数とみなし、データ $X$ との同時分布を考える。離散型の場合も含めて「密度」と書き、離散型では積分を和に読み替える。

定義 7.1（事前分布・事後分布, prior / posterior distribution）$\theta$ の密度 $\pi(\theta)$ を事前分布、$\theta$ が与えられたときの $X$ の条件付き密度 $f(x \mid \theta)$ を統計モデルとし、$(\theta, X)$ の同時密度を $\pi(\theta)f(x \mid \theta)$ とする。$X$ の周辺密度

$$m(x) = \int_\Theta f(x \mid \theta)\pi(\theta)\ d\theta$$

を周辺尤度 (marginal likelihood) という。$m(x) > 0$ となる $x$ について、$X = x$ が与えられたときの $\theta$ の条件付き密度

$$\pi(\theta \mid x) = \frac{f(x \mid \theta)\pi(\theta)}{m(x)}$$

を事後分布という。

後の式は条件付き密度の定義（第1章 定義 1.11）そのものであり、密度についてのベイズの定理である。$x$ を固定すると $f(x \mid \theta)$ は尤度 $L(\theta)$ で、$m(x)$ は $\theta$ によらないから、$\pi(\theta \mid x) \propto L(\theta)\pi(\theta)$（事後 $\propto$ 尤度 $\times$ 事前）である。ここで $g \propto h$ は、$\theta$ によらない正の定数 $c$ について $g = ch$ となることを表す。右辺を $\theta$ の関数とみて、積分が $1$ になるよう定数倍したものが事後密度である。

例 7.2（購入率）20 人中の購入者数を $X$ とし、$X \mid p \sim B(20, p)$、事前分布を一様分布 $U(0, 1)$ とする。$x = 3$ なら、$0 < p < 1$ で

これは $\operatorname{Beta}(4, 18)$ の密度の定数倍なので、事後分布は $\operatorname{Beta}(4, 18)$ である。事後平均は $4/22 = 0.182$ で、$P(p > 0.1 \mid X = 3) = 0.848$ である（数値は計算機による）。

注意 7.3（逐次更新）$X_1, X_2$ が $\theta$ を与えたもとで条件付き独立なら、$\pi(\theta \mid x_1, x_2) \propto f(x_2 \mid \theta)f(x_1 \mid \theta)\pi(\theta) \propto f(x_2 \mid \theta)\pi(\theta \mid x_1)$ である。$x_1$ による事後分布を新しい事前分布として $x_2$ で更新しても、まとめて更新しても結果は同じである。データが毎日届くなら、昨日の事後分布が今日の事前分布になる。

7.2 共役事前分布

例 7.2 では、ベータ分布（一様分布は $\operatorname{Beta}(1, 1)$）の事前分布から、ベータ分布の事後分布が得られた。

定義 7.4（共役事前分布, conjugate prior）統計モデル $f(x \mid \theta)$ に対し、事前分布の族 $\mathcal{F}$ に属するどの事前分布とどの観測値についても事後分布が $\mathcal{F}$ に属するとき、$\mathcal{F}$ を共役事前分布族という。

定理 7.5（共役事前分布）$a, b, \alpha, \beta, \tau_0, \sigma > 0$、$\mu_0 \in \mathbb{R}$ とする。ガンマ分布 $\operatorname{Gamma}(\alpha, \beta)$ の第 2 パラメータは、第1章と同じく率である。

1. （ベータ–二項）$X \mid p \sim B(n, p)$、$p \sim \operatorname{Beta}(a, b)$ ならば、$p \mid X = x \sim \operatorname{Beta}(a + x, b + n - x)$。
2. （ガンマ–ポアソン）$X_1, \dots, X_n$ が $\lambda$ を与えたもとで独立に $\operatorname{Po}(\lambda)$ に従い、$\lambda \sim \operatorname{Gamma}(\alpha, \beta)$ ならば、$\lambda \mid x \sim \operatorname{Gamma}\left(\alpha + \sum_{i} x_i, \beta + n\right)$。
3. （正規–正規、分散既知）$X_1, \dots, X_n$ が $\mu$ を与えたもとで独立に $N(\mu, \sigma^2)$ に従い（$\sigma^2$ は既知）、$\mu \sim N(\mu_0, \tau_0^2)$ ならば、$\mu \mid x \sim N(\mu_n, \tau_n^2)$。ここで

$$\frac{1}{\tau_n^2} = \frac{1}{\tau_0^2} + \frac{n}{\sigma^2}, \qquad \mu_n = \tau_n^2\left(\frac{\mu_0}{\tau_0^2} + \frac{n\bar{x}}{\sigma^2}\right)$$

証明. いずれも、尤度 $\times$ 事前密度を計算し、既知の分布の密度の定数倍になることを確かめる。密度の定数倍で積分が $1$ になるものはその密度自身に限るので、事後分布はその分布である。

1. $0 < p < 1$ で $\pi(p \mid x) \propto p^{x}(1 - p)^{n-x} \cdot p^{a-1}(1 - p)^{b-1} = p^{a+x-1}(1 - p)^{b+n-x-1}$。

2. $\lambda > 0$ で $\pi(\lambda \mid x) \propto \prod_{i} e^{-\lambda}\lambda^{x_i} \cdot \lambda^{\alpha-1}e^{-\beta\lambda} = \lambda^{\alpha + \sum_i x_i - 1}e^{-(\beta + n)\lambda}$（$1/x_i!$ は $\lambda$ によらない）。

3. $\sum_{i}(x_i - \mu)^2 = n(\mu - \bar{x})^2 + \sum_{i}(x_i - \bar{x})^2$ の第 2 項は $\mu$ によらないから

指数の中の $\mu^2$ の係数は $-\frac{1}{2}\left(\frac{n}{\sigma^2} + \frac{1}{\tau_0^2}\right) = -\frac{1}{2\tau_n^2}$、$\mu$ の係数は $\frac{n\bar{x}}{\sigma^2} + \frac{\mu_0}{\tau_0^2} = \frac{\mu_n}{\tau_n^2}$ なので、指数の中は $-\frac{(\mu - \mu_n)^2}{2\tau_n^2}$ に $\mu$ によらない定数を加えたものである。よって事後分布は $N(\mu_n, \tau_n^2)$ である。$\square$

事後平均を書き直すと、どれも事前平均と最尤推定値の重み付き平均になる：

事前分布 $\operatorname{Beta}(a, b)$ は「$a + b$ 人分のデータ（購入 $a$ 人）」、$\operatorname{Gamma}(\alpha, \beta)$ は「$\beta$ 期間に $\alpha$ 件」と同じ重みをもつ。正規–正規では、分散の逆数（精度）が足し算になる：事後の精度 $=$ 事前の精度 $+$ データの精度。$n \to \infty$ ではデータの重みが $1$ に近づき、事前分布の影響は消える。

例 7.6（事前情報を入れる）例 7.2 で、過去のページの実績から事前分布を $\operatorname{Beta}(2, 18)$（平均 $0.1$、20 人分の重み）とすると、事後分布は $\operatorname{Beta}(5, 35)$ で、事後平均は $0.5 \cdot 0.1 + 0.5 \cdot 0.15 = 0.125$ である。

例 7.7（1 日の注文数）ある商品の 1 日の注文数を $\operatorname{Po}(\lambda)$ とし、類似商品の実績から $\lambda \sim \operatorname{Gamma}(4, 2)$（平均 $2$、分散 $1$）とする。5 日間の注文数が $3, 1, 4, 2, 5$（合計 $15$）なら、事後分布は $\operatorname{Gamma}(19, 7)$ で、事後平均は $19/7 = \frac{2}{7} \cdot 2 + \frac{5}{7} \cdot 3 = 2.71$、事後標準偏差は $\sqrt{19}/7 = 0.62$ である。

例 7.8（A/B テストの効果の縮小）ある施策による売上の改善率の推定値が $\bar{x} = 2.0$（%）、標準誤差が $1.0$ だった（$z = 2.0$、両側 $p$ 値 $0.046$）。正規近似で $\bar{x} \mid \mu \sim N(\mu, 1.0^2)$ とみなし、過去の多数の実験での効果の分布から事前分布を $\mu \sim N(0, 0.5^2)$ とする。定理 7.5 の 3（$n = 1$）より、事後の精度は $1/0.25 + 1/1 = 5$ で、事後分布は $N(0.4, 0.2)$ である。推定値は 5 分の 1 に縮小され、$P(\mu > 0 \mid \bar{x}) = \Phi(0.4/\sqrt{0.2}) = 0.81$ となる。

7.3 ベイズ推定量

事後分布から一つの値を報告するとき、どの値を選ぶべきかは、外れたときの損失で決まる。

定義 7.9（ベイズ推定量, Bayes estimator）推定値 $a$ を報告し、真の値が $\theta$ だったときの損失を表す関数 $L(\theta, a) \geq 0$ を損失関数 (loss function) という。観測値 $x$ のもとでの事後期待損失

$$\rho(a \mid x) = E[L(\theta, a) \mid X = x] = \int_\Theta L(\theta, a)\pi(\theta \mid x)\ d\theta$$

を最小にする $a$ を $\delta(x)$ とするとき、$\delta$ をベイズ推定量という。

定理 7.10 $\theta$ は 1 次元とし、$x$ を固定する。

1. （二乗損失）$L(\theta, a) = (\theta - a)^2$、$E[\theta^2 \mid X = x] < \infty$ ならば、$\rho(a \mid x) = \operatorname{Var}(\theta \mid X = x) + (a - E[\theta \mid X = x])^2$ であり、ベイズ推定量は事後平均 $E[\theta \mid X = x]$ である（最小点は一意）。
2. （絶対値損失）$L(\theta, a) = \lvert \theta - a \rvert$、$E[\lvert \theta \rvert \mid X = x] < \infty$ ならば、事後分布の中央値 $m$（$P(\theta \leq m \mid X = x) \geq 1/2$ かつ $P(\theta \geq m \mid X = x) \geq 1/2$ をみたす数）は $\rho(a \mid x)$ を最小にする。

証明. 確率と期待値はすべて事後分布についてのものとし、「$\mid X = x$ 」を省く。

1. $\bar{\theta} = E[\theta]$ とおくと $(\theta - a)^2 = (\theta - \bar{\theta})^2 + 2(\theta - \bar{\theta})(\bar{\theta} - a) + (\bar{\theta} - a)^2$ で、中央の項の期待値は $0$ だから $\rho(a) = \operatorname{Var}(\theta) + (a - \bar{\theta})^2$。これは $a = \bar{\theta}$ でだけ最小になる。

2. $a > m$ とする。$\theta \leq m$ なら $\lvert \theta - a \rvert - \lvert \theta - m \rvert = a - m$、$\theta > m$ なら三角不等式より $\lvert \theta - a \rvert - \lvert \theta - m \rvert \geq -(a - m)$ なので

$a < m$ なら、$\theta \geq m$ で差は $m - a$、$\theta < m$ で差は $-(m - a)$ 以上なので、$\rho(a) - \rho(m) \geq (m - a)\bigl(2P(\theta \geq m) - 1\bigr) \geq 0$。$\square$

注意 7.11（ベイズリスク）推定量 $\delta$ の頻度論的なリスク $R(\theta, \delta) = E[L(\theta, \delta(X)) \mid \theta]$（二乗損失なら第3章の平均二乗誤差）を事前分布で平均したベイズリスク $\int_\Theta R(\theta, \delta)\pi(\theta)\ d\theta$ は、積分の順序を交換すると $E[\rho(\delta(X) \mid X)]$ に等しい。よってベイズ推定量はベイズリスクも最小にする。二乗損失では、これは「$X$ の関数で $\theta$ を最もよく近似するのは $E[\theta \mid X]$ である」という事実（11 確率論 第5章 定理 5.3）にほかならない。

例 7.12 例 7.2 の事後分布 $\operatorname{Beta}(4, 18)$ では、事後平均は $0.182$、事後中央値は $0.172$、事後密度が最大になる点（事後最頻値, MAP 推定値）は $3/20 = 0.15$ である。一様事前分布では事後密度が尤度に比例するので、MAP 推定値は最尤推定値と一致する。事後平均 $(x + 1)/(n + 2)$ は不偏ではないが、$p$ が $1/2$ に近ければ最尤推定量より平均二乗誤差が小さい（問題 7.4）。

7.4 信用区間と信頼区間

定義 7.13（信用区間, credible interval）$0 < \alpha < 1$ とする。観測値 $x$ から決まる区間 $C(x)$ で $P(\theta \in C(x) \mid X = x) = 1 - \alpha$ となるものを、$\theta$ の $100(1 - \alpha)$ % 信用区間という。事後分布の下側 $\alpha/2$ 点と $1 - \alpha/2$ 点を両端とするものを等裾信用区間、事後密度がある値以上となる $\theta$ 全体として作るもの（事後密度が単峰なら区間になり、同じ確率の区間のうちで最も短い）を最高事後密度区間（HPD 区間）という。

第4章の信頼区間とは、確率の意味が違う。

例 7.14（購入率の区間）20 人中 3 人（例 7.2）について計算すると次のようになる（数値は計算機による）。

• 一様事前分布による 95% 等裾信用区間は $[0.054, 0.363]$、HPD 区間は $[0.041, 0.340]$。
• 正規近似による 95% 信頼区間（ワルド区間）$\hat{p} \pm 1.96\sqrt{\hat{p}(1 - \hat{p})/n}$ は $[-0.006, 0.306]$。下端が負であり、小標本では使えない。
• 二項分布の確率を正確に使う 95% 信頼区間（クロッパー–ピアソン区間。両端はベータ分布の分位点で表せる）は $[0.032, 0.379]$。

注意 7.15（両者が近づく場合）定理 7.5 の 3 で $\tau_0 \to \infty$ とすると、事後分布は $N(\bar{x}, \sigma^2/n)$ に近づく。これは $\mathbb{R}$ 上の「一様分布」$\pi(\mu) = 1$ を形式的に使った事後分布である（積分が無限大なので確率分布ではなく、非正則な事前分布 (improper prior) という。使うときは尤度 $\times$ 事前の積分が有限であることを確かめる）。その 95% 等裾信用区間 $\bar{x} \pm 1.96\sigma/\sqrt{n}$ は、分散既知の正規平均の 95% 信頼区間と数値として一致する。より一般に、データが i.i.d. で統計モデルが正しく（真の分布がモデルに含まれ）、パラメータが有限次元で正則条件がみたされ、真の値がパラメータ空間の内部にあり、事前密度が真の値の近くで連続かつ正なら、$n \to \infty$ で事後分布は最尤推定量を中心とする正規分布（分散はフィッシャー情報量の逆数を $n$ で割ったもの）に近づき、信用区間は近似的に信頼区間としても使える（ベルンシュタイン–フォン・ミーゼスの定理。主張のみ。Gelman et al. の Bayesian Data Analysis を参照）。小標本や、パラメータが境界に近い場合（購入率が $0$ に近いなど）には両者は一致しない。モデルが誤っている場合も、事後分布の広がりは最尤推定量の実際のばらつきと一般には一致せず、信用区間を信頼区間として読むことはできない。

7.5 事後予測分布

次の 100 人のうち何人が購入するかを予測したい。最尤推定値を代入した $B(100, 0.15)$ で予測すると、$p$ の不確かさを無視することになる。

定義 7.16（事後予測分布, posterior predictive distribution）将来の観測 $\tilde{X}$ が、$\theta$ を与えたもとで $X$ と条件付き独立で、条件付き密度 $g(\tilde{x} \mid \theta)$ をもつとする。$X = x$ が与えられたときの $\tilde{X}$ の条件付き密度

$$p(\tilde{x} \mid x) = \int_\Theta g(\tilde{x} \mid \theta)\pi(\theta \mid x)\ d\theta$$

を事後予測分布という。

実際、$(\theta, X, \tilde{X})$ の同時密度は $\pi(\theta)f(x \mid \theta)g(\tilde{x} \mid \theta)$ なので、$\theta$ について積分して $m(x)$ で割ればこの式を得る。

命題 7.17（ベータ二項分布）事後分布が $\operatorname{Beta}(a', b')$ で、$\tilde{X} \mid p \sim B(m, p)$ ならば、$k = 0, 1, \dots, m$ について

$$P(\tilde{X} = k \mid X = x) = \binom{m}{k}\frac{B(a' + k, b' + m - k)}{B(a', b')}$$

であり（ベータ二項分布）、$s = a' + b'$ とおくと $E[\tilde{X} \mid X = x] = ma'/s$、$\operatorname{Var}(\tilde{X} \mid X = x) = \dfrac{ma'b'(s + m)}{s^2(s + 1)}$ である。

証明. 確率関数は

である。以下すべて $X = x$ のもとで考え、$p$ を与えたときの $\tilde{X}$ の条件付き平均 $mp$ と条件付き分散 $mp(1 - p)$ に全期待値・全分散の公式（第1章 命題 1.12）を使う。ベータ分布の平均と分散（第1章 1.3 節の表）より $E[p] = a'/s$、$\operatorname{Var}(p) = a'b'/(s^2(s + 1))$ で、$E[p(1 - p)] = E[p] - E[p]^2 - \operatorname{Var}(p) = a'b'/(s(s + 1))$ だから

分散の第 1 項は $p$ がわかっていても残る揺らぎ、第 2 項は $p$ の不確かさによる揺らぎである。

例 7.18（次の 100 人）事後分布 $\operatorname{Beta}(4, 18)$（例 7.2）のもとで、次の $m = 100$ 人の購入者数 $\tilde{X}$ の事後予測分布は平均 $18.2$、標準偏差 $8.9$ で、$P(\tilde{X} \leq 5 \mid x) = 0.048$、$P(\tilde{X} \geq 35 \mid x) = 0.050$ より、90% 以上の確率で $6$ 人以上 $34$ 人以下である。最尤推定値を代入した $B(100, 0.15)$ では標準偏差 $3.6$、同じ作り方（両側の確率がそれぞれ 5% 以下）の範囲は $9$ 人から $21$ 人で、ずっと狭い。代入による予測は、在庫や人員の計画で需要の揺らぎを過小評価させる。なお $m = 1$ とすると次の 1 人が購入する確率は $a'/(a' + b')$ であり、一様事前分布ではこれは $(x + 1)/(n + 2)$ となる（ラプラスの継起則）。

7.6 マルコフ連鎖モンテカルロ法

共役でないモデルでは、$m(x)$ が $k$ 次元の積分になり、$k$ が大きいと数値積分も難しい。しかし事後分布に従う乱数 $\theta_1, \dots, \theta_N$ が得られれば、事後平均は標本平均で、信用区間は標本分位点で近似できる。密度の比だけを使って、事後分布を定常分布にもつマルコフ連鎖を作り、長く走らせるのがマルコフ連鎖モンテカルロ法 (Markov chain Monte Carlo, MCMC) である。この節では、目標とする密度（ベイズ統計では事後密度）を $\pi$ と書き、$\mathbb{R}^k$ 上の関数とみなす（$\Theta$ の外では $0$）。$\pi$ は定数倍を除いてわかればよい。11 確率論 第6章 の可算な状態空間の推移確率の代わりに、推移核を使う。

定義 7.19（推移核・定常分布・詳細つり合い）$\theta \in \mathbb{R}^k$ と集合 $A \subset \mathbb{R}^k$ に確率 $K(\theta, A)$ を対応させ、各 $\theta$ について $K(\theta, \cdot)$ が確率分布であるものを推移核という。$\theta_{t+1}$ の条件付き分布が、過去によらず $K(\theta_t, \cdot)$ で与えられる列 $(\theta_t)$ をマルコフ連鎖という。確率密度 $\pi$ がすべての $A$ について

$$\int \pi(\theta)K(\theta, A)\ d\theta = \int_A \pi(\theta)\ d\theta$$

をみたすとき、$\pi$ を $K$ の定常分布という。すべての $A, B$ について

$$\int_A \pi(\theta)K(\theta, B)\ d\theta = \int_B \pi(\theta)K(\theta, A)\ d\theta$$

となるとき、$K$ は $\pi$ について詳細つり合い (detailed balance) をみたすという。

定常分布の式は「$\theta_t \sim \pi$ ならば $\theta_{t+1} \sim \pi$」を意味する。詳細つり合いの式で $B = \mathbb{R}^k$ とおくと、$K(\theta, \mathbb{R}^k) = 1$ より定常分布の式を得る。つまり、詳細つり合いをみたす $\pi$ は定常分布である（可算な場合は 11 確率論 第6章 定義 6.22 の後の議論）。

定義 7.20（メトロポリス–ヘイスティングス法, Metropolis–Hastings algorithm）各 $\theta$ について $q(\theta, \cdot)$ が確率密度であるもの（提案分布）をとり、$\pi(\theta)q(\theta, \theta') > 0$ のとき

$$\alpha(\theta, \theta') = \min\left\lbrace 1, \frac{\pi(\theta')q(\theta', \theta)}{\pi(\theta)q(\theta, \theta')} \right\rbrace$$

とおく（それ以外では $\alpha = 1$ とする）。$\pi(\theta_0) > 0$ となる $\theta_0$ から始め、現在の値 $\theta_t$ から次の手順で $\theta_{t+1}$ を作る。(i) $\theta' \sim q(\theta_t, \cdot)$ を生成する。(ii) 確率 $\alpha(\theta_t, \theta')$ で受理して $\theta_{t+1} = \theta'$ とし、受理しなければ（棄却）$\theta_{t+1} = \theta_t$ とする。

$\alpha$ には $\pi$ の比しか現れないので、正規化定数は要らない。提案が対称、すなわち $q(\theta, \theta') = q(\theta', \theta)$ なら（$\theta' = \theta + \varepsilon$、$\varepsilon \sim N(0, s^2I)$ など）$\alpha = \min\lbrace 1, \pi(\theta')/\pi(\theta) \rbrace$ となる。これが最初のメトロポリス法（1953 年）で、ヘイスティングス（1970 年）が非対称な提案に拡張した。この手順の推移核は

である。第 1 項は提案が受理されて $A$ に移る確率、$r(\theta)$ は棄却されてその場に留まる確率である。

定理 7.21（メトロポリス–ヘイスティングス法の正しさ）定義 7.20 の推移核 $K$ は $\pi$ について詳細つり合いをみたす。特に、$\pi$ は $K$ の定常分布である。

証明. $g(\theta, \theta') = \pi(\theta)q(\theta, \theta')\alpha(\theta, \theta')$ とおく。$\pi(\theta)q(\theta, \theta') > 0$ なら

であり、$\pi(\theta)q(\theta, \theta') = 0$ なら両辺とも $0$ なので、この式はつねに成り立つ。右辺は $\theta$ と $\theta'$ について対称だから $g(\theta, \theta') = g(\theta', \theta)$ である。$K$ の式より

右辺の第 1 項は、積分の順序を交換し（被積分関数が非負なので交換できる。厳密には 06 測度と積分 第5章 定理 5.6（トネリの定理））、$g$ の対称性を使うと $\int_B\left(\int_A g(\theta, \theta')\ d\theta'\right)d\theta$ に等しい。第 2 項は $A$ と $B$ について対称である。よって左辺は $A$ と $B$ を入れ替えても変わらない。後半は定義 7.19 の後の注意による。$\square$

注意 7.22（定常分布にもつことと、収束すること）定理 7.21 が保証するのは「$\theta_t \sim \pi$ ならその後もずっと $\pi$ に従う」ことだけである。任意の初期値から $\theta_t$ の分布が $\pi$ に近づくことや、時間平均 $\frac{1}{N}\sum_{t=1}^{N}h(\theta_t)$ が $\int h\pi$ に収束することには、別の条件が要る。状態空間が可算なら、連鎖が既約（どの状態からどの状態へも正の確率で到達できる）であれば、定常分布 $\pi$ をもつことから連鎖は正再帰的で、定常分布は $\pi$ だけである（11 確率論 第6章 定理 6.17）。さらに $\sum \lvert h \rvert\pi < \infty$ なら、任意の初期値から時間平均は確率 1 で $\sum h\pi$ に収束し（同 定理 6.21）、連鎖が非周期的なら任意の初期値から分布が $\pi$ に収束する（11 確率論 第6章 定理 6.20）。$\mathbb{R}^k$ 上でも、既約性・非周期性にあたる条件のもとで同様の定理が成り立ち、たとえば $\pi$ が $\mathbb{R}^k$ 上で連続かつ正で、提案が正規分布のランダムウォークなら条件はみたされる（一般の状態空間のマルコフ連鎖の理論による。主張のみ。Gelman et al. を参照）。条件が崩れる例：目標 $\pi$ が $[-3, -2] \cup [2, 3]$ 上の一様分布、提案が $\theta' \sim U(\theta - 1, \theta + 1)$ なら、$-2.5$ から出発した連鎖は、$[-3, -2]$ の外への提案が $\pi(\theta') = 0$ のためすべて棄却されるので、$[-3, -2]$ から出られない。$\pi$ は定常分布だが（$[-3, -2]$ 上の一様分布も定常分布である）、連鎖の分布は $\pi$ に近づかない。

例 7.23（共役でない事前分布）例 7.2 で、購入率 $p$ の代わりに対数オッズ $\theta = \log(p/(1 - p))$ に事前分布 $N(-2, 1)$ をおく（$p$ の事前中央値は $1/(1 + e^{2}) = 0.12$）。事後密度は

で、共役ではない。提案 $\theta' = \theta + \varepsilon$（$\varepsilon \sim N(0, 1)$）のメトロポリス法で乱数を作る。桁あふれを避けるため対数で計算し、一様乱数 $U$ について $\log U < \log \pi(\theta') - \log \pi(\theta)$ のとき受理する（確率 $\min\lbrace 1, \pi(\theta')/\pi(\theta) \rbrace$ での受理と同じである）。

import numpy as np

def log_post(t, x=3, n=20):      # θ = log(p/(1-p)) の対数事後密度（定数を除く）
    return x * t - n * np.log1p(np.exp(t)) - (t + 2) ** 2 / 2

rng = np.random.default_rng(0)
t, chain, acc = 0.0, [], 0
for k in range(50000):
    s = t + rng.normal(0, 1.0)   # 対称な提案（ランダムウォーク）
    if np.log(rng.uniform()) < log_post(s) - log_post(t):
        t, acc = s, acc + 1      # 受理（棄却なら t に留まる）
    chain.append(t)
p = 1 / (1 + np.exp(-np.array(chain[5000:])))   # 最初の 5000 回を捨て、p に戻す
print(f"受理率 {acc / 50000:.3f}")
print(f"事後平均 {p.mean():.4f}")
print(f"95% 信用区間 [{np.quantile(p, 0.025):.4f}, {np.quantile(p, 0.975):.4f}]")

受理率 0.526
事後平均 0.1442
95% 信用区間 [0.0462, 0.2965]

格子上の数値積分で求めた値（事後平均 $0.1442$、95% 等裾信用区間 $[0.0458, 0.2954]$）とよく一致する。最初の部分を捨てる（バーンイン, burn-in）のは、初期値 $0$ の影響が残る部分を除くためである。連続する値は互いに相関しているので、45000 個の値は 45000 個の独立な標本ほどの精度をもたない。

パラメータが多いモデルでは、ほかの成分を固定したときの 1 成分の条件付き分布（完全条件付き分布, full conditional distribution）が簡単なことが多い。$\theta = (\theta_1, \dots, \theta_k)$ から第 $j$ 成分を除いたものを $\theta_{-j}$ と書く。

定義 7.24（ギブスサンプリング, Gibbs sampling）現在の値 $\theta$ の第 $j$ 成分を、完全条件付き分布 $\pi(\theta_j \mid \theta_{-j})$ から生成した値で置き換える操作を、$j = 1, \dots, k$ の順に行うことを 1 回の更新とする。

定理 7.25 各成分の置き換えは $\pi$ を不変にする。すなわち $\theta \sim \pi$ ならば、第 $j$ 成分を置き換えた値も $\pi$ に従う。したがって $\pi$ はギブスサンプリングの定常分布である。

証明. $\theta \sim \pi$ なら、$\theta_{-j}$ は周辺密度 $\pi_{-j}(\theta_{-j}) = \int \pi(\theta)\ d\theta_j$ に従う。置き換えた値 $\theta'_j$ は $\theta_{-j}$ を与えたもとで密度 $\pi(\theta'_j \mid \theta_{-j})$ に従うので、$(\theta'_j, \theta_{-j})$ の同時密度は $\pi_{-j}(\theta_{-j})\pi(\theta'_j \mid \theta_{-j}) = \pi(\theta'_j, \theta_{-j})$（条件付き密度の定義）である。$\pi$ を不変にする操作を続けて行っても $\pi$ は不変である。$\square$

成分の置き換えは、提案を $\pi(\theta'_j \mid \theta_{-j})$ としたメトロポリス–ヘイスティングス法で、形式的に受理確率を計算すると $1$ になるものとみることもできる。収束については注意 7.22 と同じ注意が要り、成分の間の相関が強いと、1 成分ずつしか動けないために連鎖の動きが遅くなる（問題 7.7）。

例 7.26（平均と分散が未知の正規モデル）$X_1, \dots, X_n$ が $(\mu, \tau)$ を与えたもとで独立に $N(\mu, 1/\tau)$ に従い（$\tau$ は精度）、事前分布は独立に $\mu \sim N(\mu_0, \tau_0^2)$、$\tau \sim \operatorname{Gamma}(\alpha_0, \beta_0)$ とする。同時事後分布はよく知られた分布ではないが、完全条件付き分布は

である。前者は定理 7.5 の 3 で $\sigma^2 = 1/\tau$ としたもの、後者は $\pi(\tau \mid \mu, x) \propto \tau^{n/2}e^{-\tau\sum_i (x_i - \mu)^2/2} \cdot \tau^{\alpha_0 - 1}e^{-\beta_0\tau}$ による。この 2 つを交互に生成すればよい。

7.7 階層モデルと縮小推定（紹介）

店舗・広告・地域など $J$ 個のグループのパラメータを同時に推定したい。グループ $j$ の推定値 $y_j$ が $y_j \mid \theta_j \sim N(\theta_j, \sigma_j^2)$（$\sigma_j$ は既知の標準誤差）に従うとする。グループごとに別々に推定する（$\hat{\theta}_j = y_j$）と、データの少ないグループの推定値は大きくばらつく。全体を一つにまとめると、グループの本当の違いが消える。階層モデル (hierarchical model) ではその中間をとり、$\theta_j$ 自体が共通の分布 $N(\mu, \eta^2)$ から独立に生じると考える。$\mu, \eta$ を既知とすれば、定理 7.5 の 3 より