統計学 学習ノート

はじめに

「事前分布なんて人が決めるものなのに、それで『最適』と言えるの？」

準1級でベイズ法に触れたとき、こう感じた人もいるかもしれません。事前分布を変えれば、答えも変わります。それでは、理論としていい加減ではないか、という疑問です。

例えば、天気予報を思い浮かべてください。降水確率の予報は、人によって少しずつ違います。それでも、「傘を持つかどうか」という判断には、はっきりした考え方があります。濡れたときの困り具合と、傘を持ち歩く手間を比べて、損が最も小さくなるほうを選ぶのです。

実は、ベイズ法も同じです。事前分布と損失を決めれば、最もよい推定量が定理として導けます。さらに、事前分布を別のデータから学ぶ仕組みまで作れます。ここまでの1級の記事で身につけた、尤度、十分統計量、MCMCの知識は、すべてこの1本にまとまります。

そこで登場するのが、今回の主役 事後分布 です。その計算を支える ギブスサンプリング と、データ全体で情報を分け合う 階層ベイズモデル も、一緒に導いていきます。

前回の記事では、MCMCの正しさを証明しました。今回の記事では、事前分布 から事後分布を導く手順、ベイズ推定量の証明、ギブスサンプリングの理論、階層ベイズモデルの計算について解説します。統計数理の最後の記事です。

事前分布から事後分布を導く

はじめに、事後分布の導き方について解説します。

母数 $\theta$ の不確かさを表す分布を 事前分布 $\pi(\theta)$ と呼び、データ $x$ を見たあとの分布を 事後分布 $\pi(\theta \mid x)$ と呼びます。ベイズの定理から、次の関係がありました。

$$\pi(\theta \mid x) = \frac{f(x \mid \theta)\pi(\theta)}{\int f(x \mid \theta')\pi(\theta')\, d\theta'} \propto f(x \mid \theta)\pi(\theta)$$

1級の答案では、次の手順で書くと迷いません。

【事後分布を導く手順】

1. 尤度と事前分布を掛ける
2. $\theta$ を含まない因子をすべて捨てる
3. 残った式の形（カーネル）から、分布の名前を見抜く
4. 分布の名前から正規化定数を戻す

ジグソーパズルで、絵柄の特徴的なピースだけを見て何の絵か当てるようなものです。定数は最後に、分布の公式から自動的に戻ってきます。

例題：ポアソン分布とガンマ事前分布

$x_1, \dots, x_n$ が独立に平均 $\lambda$ のポアソン分布に従い、事前分布は密度が $\lambda^{a-1}e^{-b\lambda}$ に比例するガンマ分布 $\text{Ga}(a, b)$ とする。$s = \sum x_i$ とする。〔1〕事後分布を導け。〔2〕事後平均が事前平均と標本平均の重みつき平均であることを示せ。〔3〕次の1個の観測 $X$ の予測分布を導け。〔4〕$a = 3$、$b = 1$、$n = 5$、$s = 20$ のとき、事後平均、MAP推定値、$X$ の予測平均と予測分散を求めよ。

【解答】

〔1〕尤度の $\lambda$ を含む部分は $\lambda^{s}e^{-n\lambda}$ である。事前分布を掛けて、次を得る。

$$\pi(\lambda \mid x) \propto \lambda^{s}e^{-n\lambda} \cdot \lambda^{a-1}e^{-b\lambda} = \lambda^{a+s-1}e^{-(b+n)\lambda}$$

これは $\text{Ga}(a + s, b + n)$ のカーネルである。

〔2〕事後平均は、次のように変形できる。

$$\frac{a + s}{b + n} = \frac{b}{b + n} \cdot \frac{a}{b} + \frac{n}{b + n} \cdot \bar{x}$$

$a/b$ は事前平均なので、重みが $b$ と $n$ の重みつき平均である。

〔3〕$a' = a + s$、$b' = b + n$ とおく。$\lambda$ について積分すると、ガンマ関数の積分から次を得る。

$$P(X = k \mid x) = \int_0^{\infty} \frac{e^{-\lambda}\lambda^{k}}{k!} \cdot \frac{b'^{a'}}{\Gamma(a')}\lambda^{a'-1}e^{-b'\lambda}\, d\lambda = \frac{\Gamma(a' + k)}{k!\,\Gamma(a')}\left(\frac{b'}{b' + 1}\right)^{a'}\left(\frac{1}{b' + 1}\right)^{k}$$

これは負の二項分布である。

〔4〕事後分布は $\text{Ga}(23, 6)$ である。事後平均は $23/6 = 3.83$、MAP推定値は $22/6 = 3.67$ である。予測平均は $E(\lambda \mid x) = 3.83$ である。予測分散は、全分散の公式から次のようになる。

$$V(X \mid x) = E(\lambda \mid x) + V(\lambda \mid x) = \frac{23}{6} + \frac{23}{36} = 4.47$$

予測分散が平均より大きいのは、$\lambda$ の不確かさの分が加わるからです。ここが、$\hat{\lambda}$ を代入しただけの予測との違いです。

もう1つの基本形は、正規分布の平均です。$\sigma^2$ を既知として、$\mu$ の事前分布を $N(m_0, \tau_0^2)$ とします。指数の中の $\mu$ の2次式を、平方完成します。

$$-\frac{1}{2}\left\{\frac{n(\mu - \bar{y})^2}{\sigma^2} + \frac{(\mu - m_0)^2}{\tau_0^2}\right\} = -\frac{1}{2}\left(\frac{n}{\sigma^2} + \frac{1}{\tau_0^2}\right)(\mu - m_1)^2 + \text{定数}$$

ここから、事後の精度は2つの精度の和で、事後平均 $m_1$ は精度で重みをつけた平均だと分かります。準1級の公式が、平方完成の1行から出ました。

次に、一般的な性質を2つ確認します。

【事後分布の一般的な性質】

1. 事後分布は十分統計量だけに依存する： 分解定理 $f(x \mid \theta) = g(T(x), \theta)h(x)$ の $h(x)$ が約分で消える
2. データを分けて更新しても結果は同じ： $x_1$ の事後分布を事前分布として $x_2$ で更新すると、$f(x_2 \mid \theta)f(x_1 \mid \theta)\pi(\theta)$ に比例する

最後に、事前分布の選び方です。情報がないときの候補に、ジェフリーズの事前分布 があります。フィッシャー情報量 $I(\theta)$ を使って、次のように定めます。

$$\pi(\theta) \propto \sqrt{I(\theta)}$$

利点は、母数の変換に左右されないことです。$\phi = g(\theta)$ と変換すると、$I(\phi) = I(\theta)(d\theta/d\phi)^2$ です。したがって、$\sqrt{I(\phi)} = \sqrt{I(\theta)}\,|d\theta/d\phi|$ となり、これは密度の変数変換の公式そのものです。どちらの母数で「無情報」と言っても同じ結論になるのです。

例えば、二項分布なら $I(p) = n/\{p(1-p)\}$ なので、事前分布は $p^{-1/2}(1-p)^{-1/2}$ に比例します。これは $\text{Beta}(1/2, 1/2)$ です。10回中7回成功なら、事後分布は $\text{Beta}(7.5, 3.5)$ で、事後平均は0.682です。

補足（ここは読み飛ばしてOK）： ポアソン分布なら $\lambda^{-1/2}$ となり、積分が発散します。こうした事前分布を非正則と呼びます。非正則でも、事後分布の積分が有限なら使えます。答案では、事後分布が確率分布になることを確かめます。

試験のポイント：答案では「$\theta$ を含む部分だけに注目する」と宣言してから比例の記号を使います。$a$ や $b$ などの定数を、最後に分布の名前で回収します。

以上が、事後分布の導き方についての解説です。

ベイズ推定量：損失から最適な推定量を導く

次に、ベイズ推定量について説明します。

推定値 $a$ を報告したときの痛手を、損失関数 $L(\theta, a)$ で表します。ベイズ推定量 は、事後分布による損失の期待値を最小にする $a$ です。

$$\delta(x) = \arg\min_{a} E\{L(\theta, a) \mid x\}$$

冒頭の傘の判断と同じ形です。この $\delta$ は、$\theta$ と $x$ の両方で平均をとったベイズリスクも最小にします。ベイズリスクは、積分の順序を入れ替えると、「$x$ ごとの事後期待損失」を $x$ の周辺分布で平均したものだからです。各 $x$ で最小にすれば、全体も最小になります。

損失の形で、答えが変わります。1つずつ導きましょう。

【損失関数とベイズ推定量】

1. 二乗損失 $(\theta - a)^2$： 事後平均
2. 絶対損失 $|\theta - a|$： 事後中央値
3. 0-1損失（ごく小さな幅で一致すれば損失0）： 事後モード（MAP）

1番の証明です。$m = E(\theta \mid x)$ とすると、次の分解ができます。

$$E\{(\theta - a)^2 \mid x\} = V(\theta \mid x) + (m - a)^2$$

第1項は $a$ を含まないので、$a = m$ で最小です。

2番の証明です。$E(|\theta - a| \mid x)$ を $a$ で微分すると、$P(\theta < a \mid x) - P(\theta > a \mid x)$ になります。これが0になるのは、両側の確率が等しいとき、つまり $a$ が事後中央値のときです。的当てに例えると、大きく外すほど罰が増える二乗損失は的の重心を、外した距離に比例する絶対損失は的の中央値を狙う、というイメージです。

例題：二項分布のベイズ推定量のリスク

$X \sim B(n, \theta)$ とし、事前分布を $\text{Beta}(a, b)$ とする。〔1〕二乗損失のベイズ推定量 $\delta$ を求めよ。〔2〕$\delta$ の平均二乗誤差を求めよ。〔3〕$a = b = \sqrt{n}/2$ のとき、平均二乗誤差が $\theta$ によらないことを示せ。〔4〕$n = 16$ のとき、この $\delta$ と最尤推定量 $X/n$ の平均二乗誤差を比べよ。

【解答】

〔1〕事後分布は $\text{Beta}(a + x, b + n - x)$ だから、$\delta = (X + a)/(n + a + b)$ である。

〔2〕分散と偏りの2乗の和である。

$$\text{MSE}(\theta) = \frac{n\theta(1 - \theta) + \{a - (a + b)\theta\}^2}{(n + a + b)^2}$$

〔3〕$a = b = \sqrt{n}/2$ を代入すると、$\{a - (a + b)\theta\}^2 = n(1/2 - \theta)^2$ だから、分子は次のようになる。

$$n\theta - n\theta^2 + \frac{n}{4} - n\theta + n\theta^2 = \frac{n}{4}$$

よって、$\text{MSE} = \dfrac{n/4}{(n + \sqrt{n})^2} = \dfrac{1}{4(\sqrt{n} + 1)^2}$ で、$\theta$ によらない。

〔4〕$n = 16$ なら $\delta = (X + 2)/20$ で、$\text{MSE} = 0.01$ である。最尤推定量の平均二乗誤差は $\theta(1 - \theta)/16$ で、$\theta = 0.5$ で0.0156、$\theta = 0.1$ で0.0056である。$0.2 < \theta < 0.8$ では $\delta$ のほうが小さい。

事前分布を置いた推定量が、最尤推定量より多くの $\theta$ で優れることがあります。逆に、0や1に近い $\theta$ では負けます。つまり、一様に最もよい推定量はありません。

補足（ここは読み飛ばしてOK）： リスクが一定のベイズ推定量は、最大リスクを最小にするミニマックス推定量になります。もし最大リスクがより小さい推定量があれば、それはベイズリスクも小さくなり、定義に矛盾するからです。

試験のポイント：どの損失でどの推定量になるかだけでなく、二乗損失と絶対損失の証明を書けるようにします。

以上が、ベイズ推定量についての説明です。

ギブスサンプリング：条件付きなら解ける

続いて、ギブスサンプリングについて解説します。

母数が複数あると、同時事後分布は共役の形になりません。しかし、他の母数を固定すれば、1つの母数の分布は知っている形になることがよくあります。この事実を使うのが ギブスサンプリング です。

母数を $\theta = (\theta_1, \dots, \theta_d)$ とします。$\theta_j$ 以外をすべて固定したときの $\theta_j$ の条件付き分布を、完全条件付き分布 と呼びます。同時事後分布の式で、$\theta_j$ を含む因子だけを残せば、これが得られます。

【ギブスサンプリングの手順】

1. 初期値を決める
2. $\theta_1$ を、他の最新の値を固定した完全条件付き分布から引く
3. $\theta_2, \dots, \theta_d$ も同様に、1つずつ順に引く
4. 2と3をくり返す

全員の席を一度に決めるのは難しくても、他の人の席を固定して1人ずつ座り直すなら簡単です。これをくり返すと、全体の席順が目標の分布に落ち着きます。

なぜ正しいのでしょうか？ 2変数で、$\theta_1$ を更新する1歩を、前回のMHアルゴリズムと見なします。提案 $q$ を完全条件付き分布 $\pi(y_1 \mid x_2)$ とします。$x = (x_1, x_2)$ から $y = (y_1, x_2)$ への採択確率の比を計算しましょう。

$$\frac{\pi(y)\, q(x \mid y)}{\pi(x)\, q(y \mid x)} = \frac{\pi(x_2)\pi(y_1 \mid x_2) \cdot \pi(x_1 \mid x_2)}{\pi(x_2)\pi(x_1 \mid x_2) \cdot \pi(y_1 \mid x_2)} = 1$$

採択確率はいつも1です。つまり、ギブスサンプリングは、提案が必ず採択されるMHアルゴリズムなのです。

具体的に見てみましょう。2値の変数 $(x, y)$ の同時確率が次の表で与えられたとします。

	$y = 0$	$y = 1$
$x = 0$	0.3	0.2
$x = 1$	0.1	0.4

条件付き確率は、$P(y = 1 \mid x = 0) = 0.4$、$P(y = 1 \mid x = 1) = 0.8$、$P(x = 1 \mid y = 0) = 1/4$、$P(x = 1 \mid y = 1) = 2/3$ です。$x$ から $y$ を引き、その $y$ から次の $x$ を引くので、$x$ の連鎖の推移確率は、途中の $y$ で場合分けして求まります。

$$P(0 \to 1) = 0.6 \times \frac{1}{4} + 0.4 \times \frac{2}{3} = \frac{5}{12}, \qquad P(1 \to 0) = 0.2 \times \frac{3}{4} + 0.8 \times \frac{1}{3} = \frac{5}{12}$$

定常分布は $(1/2, 1/2)$ で、表の周辺分布 $P(x = 0) = 0.5$ と一致しました。

せっかくなので、連続の例も見てみましょう。平均0、分散1、相関 $\rho$ の2変量正規分布では、$x \mid y \sim N(\rho y, 1 - \rho^2)$ です。$y$ を経由して $x_t$ から $x_{t+1}$ へ進むと、次のようになります。

$$x_{t+1} \mid x_t \sim N(\rho^2 x_t,\ 1 - \rho^4)$$

これは自己回帰の係数が $\rho^2$ の連鎖です。前回の公式から、分散の倍率は $(1 + \rho^2)/(1 - \rho^2)$ です。$\rho = 0.9$ なら約9.5倍で、変数どうしの相関が強いほど、連鎖の動きは遅くなります。

例題：平均と分散が未知の正規分布

データ 4, 6, 5, 7, 8 が $N(\mu, \sigma^2)$ に従う。事前分布は独立に $\mu \sim N(0, 100)$、$\sigma^2 \sim \text{IG}(1, 1)$（逆ガンマ分布、密度は $(\sigma^2)^{-a-1}e^{-b/\sigma^2}$ に比例）とする。〔1〕完全条件付き分布を導け。〔2〕$\sigma^2 = 2$ のときの $\mu$ の条件付き分布と、$\mu = 6$ のときの $\sigma^2$ の条件付き分布を求めよ。

【解答】

〔1〕同時事後分布は次の式に比例する。

$$(\sigma^2)^{-n/2}\exp\left\{-\frac{\sum(y_i - \mu)^2}{2\sigma^2}\right\}\exp\left\{-\frac{(\mu - m_0)^2}{2\tau_0^2}\right\}(\sigma^2)^{-a-1}\exp\left(-\frac{b}{\sigma^2}\right)$$

$\sigma^2$ を固定すると $\mu$ の2次式だから、正規分布である。精度は $1/\tau_0^2 + n/\sigma^2$、平均は $(m_0/\tau_0^2 + n\bar{y}/\sigma^2)$ を精度で割った値である。

$\mu$ を固定すると $\sigma^2$ の部分は $(\sigma^2)^{-(a + n/2) - 1}\exp[-\{b + \sum(y_i - \mu)^2/2\}/\sigma^2]$ だから、$\text{IG}(a + n/2,\ b + \sum(y_i - \mu)^2/2)$ である。

〔2〕$\bar{y} = 6$ である。$\sigma^2 = 2$ なら、精度は $0.01 + 2.5 = 2.51$ で、$\mu \sim N(15/2.51,\ 1/2.51) = N(5.98, 0.398)$ である。$\mu = 6$ なら、$\sum(y_i - 6)^2 = 10$ だから、$\sigma^2 \sim \text{IG}(3.5, 6)$ である。

この2つの分布を交互に20万回引くと、$\mu$ の事後平均は約5.97、95%信用区間は約 (4.4, 7.5)、$\sigma^2$ の事後平均は約3.0となりました。手で積分できない事後分布が、乱数で手に入ります。

試験のポイント：完全条件付き分布は、同時分布の式から対象の母数を含む因子だけを抜き出して求めます。ギブスサンプリングでは棄却がありません。

以上が、ギブスサンプリングについての解説です。

階層ベイズモデル：事前分布もデータから学ぶ

最後に、階層ベイズモデルを見ていきましょう。

階層ベイズモデル は、事前分布の母数にも、さらに分布を置くモデルです。冒頭の疑問への、1つの答えでもあります。事前分布を人が決める代わりに、似たグループのデータから学ぶのです。

例えば、シーズン序盤に5打数4安打の新人がいたとします。打率8割だとは、誰も思いませんよね。これまでの選手全体の打率を知っているからです。階層モデルは、この当たり前の判断を式にします。

モデルは3段で書けます。グループ $j$ の標本平均 $\bar{y}_j$ の分散を $v = \sigma^2/n$ として、次のようにします。

$$\bar{y}_j \mid \theta_j \sim N(\theta_j, v), \qquad \theta_j \mid \mu, \tau^2 \sim N(\mu, \tau^2), \qquad (\mu, \tau^2) \sim \pi(\mu, \tau^2)$$

$\mu$ と $\tau^2$ を 超母数（ハイパーパラメータ）と呼びます。$\theta_j$ は、番号を入れ替えても同じ同時分布をもつ交換可能な変数です。

$\mu$ と $\tau^2$ を固定すれば、1段目と2段目は正規分布の平均の問題そのものです。精度の重みで、$\theta_j$ の事後平均が求まります。

$$E(\theta_j \mid \bar{y}_j, \mu, \tau^2) = (1 - B)\bar{y}_j + B\mu, \qquad B = \frac{v}{v + \tau^2}$$

事後分散は $(1 - B)v$ です。$B$ を 縮小係数 と呼びます。グループの平均は、全体の平均 $\mu$ に向かって、$B$ の割合だけ縮みます。データが少なく $v$ が大きいほど、またグループの差 $\tau^2$ が小さいほど、強く縮みます。力の強い綱引きの相手に引っ張られるようなものです。

超母数は、どうやって決めるのでしょうか？ $\theta_j$ を積分で消すと、$\bar{y}_j$ の周辺分布が得られます。

$$\bar{y}_j \mid \mu, \tau^2 \sim N(\mu,\ \tau^2 + v)$$

ここがポイントです。$J$ 個のグループ平均のばらつきに、$\tau^2$ の情報が入っています。$\bar{\bar{y}}$ を $\bar{y}_j$ の平均として、$S = \sum_j (\bar{y}_j - \bar{\bar{y}})^2$ とすると、$E\{S/(J - 1)\} = \tau^2 + v$ です。この関係を使って超母数を推定値で置き換える方法を、経験ベイズ法 と呼びます。

例題：5クラスの平均の縮小推定

同じ大きさの5クラスの平均点が 58, 62, 66, 70, 74 である。各クラスの平均の分散は $v = 16$ で既知とする。〔1〕$\mu$ と $\tau^2$ をモーメント法で推定せよ。〔2〕縮小係数を求めよ。〔3〕各クラスの $\theta_j$ の推定値を求めよ。

【解答】

〔1〕$\hat{\mu} = 66$ である。$S = 64 + 16 + 0 + 16 + 64 = 160$ で、$S/(J - 1) = 40$ だから、$\hat{\tau}^2 = 40 - 16 = 24$ である。

〔2〕$\hat{B} = 16/40 = 0.4$ である。

〔3〕$\hat{\theta}_j = 66 + 0.6(\bar{y}_j - 66)$ だから、順に 61.2, 63.6, 66, 68.4, 70.8 である。

もとの平均の幅が16点あったのに、推定値の幅は9.6点に縮みました。

経験ベイズ法は、$\mu$ と $\tau^2$ の不確かさを無視しています。階層ベイズモデルでは、それも事後分布に含めて、ギブスサンプリングで計算します。$\pi(\mu) \propto 1$、$\tau^2 \sim \text{IG}(a, b)$ とすると、完全条件付き分布は次の3つです。

【階層モデルのギブスサンプリング】

1. $\theta_j \mid \mu, \tau^2, y \sim N((1 - B)\bar{y}_j + B\mu,\ (1 - B)v)$
2. $\mu \mid \theta, \tau^2 \sim N(\bar{\theta},\ \tau^2/J)$
3. $\tau^2 \mid \theta, \mu \sim \text{IG}(a + J/2,\ b + \sum_j (\theta_j - \mu)^2/2)$

どれも、これまでに導いた形の組み合わせです。

補足（ここは読み飛ばしてOK）： 3個以上の平均を同時に推定するとき、全体の二乗誤差の和の意味で、平均をそのまま使う推定量は最善ではありません。縮小するジェームズ・スタイン推定量がこれを改善します。このとき、平均へ縮める場合の係数は $(J - 3)v/S$ で、上の例なら0.2です。

試験のポイント：縮小係数は「$v$ と $\tau^2$ のどちらの分散が大きいか」で決まります。周辺分布の分散が $\tau^2 + v$ になることを、導けるようにします。

以上が、階層ベイズモデルについての解説です。

細かい計算を忘れてしまっても問題ありません。ここでは、「事前分布と損失を決めれば、あとは確率の計算だけで答えが決まる」という点だけ押さえておきましょう。

今回のまとめ

それでは、今回の内容を振り返ります。

今回はベイズ法を、導出から紹介しました。事後分布はカーネルの形を見抜いて導き、最適な推定量は損失から定理として出て、計算は完全条件付き分布の組み合わせで進みます。そして、階層ベイズモデルでは、事前分布そのものをデータから学びます。冒頭の「人が決めるのに最適と言えるのか」という疑問には、「何を前提としたかを明示して、その前提のもとで最適を証明する」と答えられます。

今回のポイントを振り返りましょう。

振り返り

• 事後分布は「尤度 × 事前分布」から、$\theta$ を含む因子だけで導く。
• ジェフリーズの事前分布は $\sqrt{I(\theta)}$ に比例し、母数の変換に左右されない。
• 二乗損失なら事後平均、絶対損失なら事後中央値がベイズ推定量である。
• ギブスサンプリングは、採択確率がいつも1のMHアルゴリズムである。
• 階層ベイズモデルでは、グループの平均が全体の平均へ $B = v/(v + \tau^2)$ の割合で縮む。

キーワードのおさらい

• 事前分布: データを見る前の母数の不確かさを表す分布
• 事後分布: データを見たあとの母数の分布。尤度 × 事前分布に比例する
• 階層ベイズモデル: 事前分布の母数にも分布を置くモデル
• ギブスサンプリング: 完全条件付き分布から1つずつ順に乱数を引くMCMC

ここまでお疲れさまでした。これで、1級の統計数理の記事はすべて終わりです。確率の公理から出発して、分布、推定、検定、回帰、分割表、順位、不完全データ、シミュレーション、ベイズ法と、理論を自分の手で導く道を歩いてきました。「公式を使う」から「公式を作る」へ、大きな一歩を踏み出したはずです。

さて、ここまでの理論は、データが手元にあることから話を始めていました。しかし、実際の仕事で最初に悩むのは、そのデータをどう集めるかです。実験をするのか、観察するのか、調査するのか。この選び方を誤ると、どんな理論でも結論を救えません。

次回からは、1級の統計応用に入ります。最初は、研究の種類として、実験・観察・調査について紹介します。

それでは、次の記事へ進んでいきましょう。