統計学 学習ノート

はじめに

母数を確率で語れば、0.5より大きい確率を素直に計算できる

「『母数が0.5より大きい確率は87%です』って、言ってはいけないの？」

2級で信頼区間を学んだとき、注意がありましたよね。95%信頼区間は、「母数がこの区間に入る確率が95%」という意味ではない、という話です。母数は決まった値で、確率的に動くのは区間のほうだからです。正しい説明ですが、少し回りくどく感じたかもしれません。

例えば、初めて入ったラーメン店がおいしかったとしましょう。「この店は当たりの可能性が高い」と考えますよね。もう1回行っておいしければ、確信はさらに強まります。私たちは日常で、知りたい値について確率で語り、データを見るたびに更新しています。

実は、この考え方をそのまま統計学にしたものがあります。母数を確率変数として扱えば、「母数が0.5より大きい確率」を素直に計算できるのです。

そこで登場するのが、今回の主役 ベイズ法 です。出発点は、準1級の最初の記事で学んだベイズの定理です。

前回の記事で、BICの導出にベイズの考え方が少しだけ出てきました。今回の記事では、事前分布から事後分布を求める方法、共役事前分布、そして計算を支えるマルコフ連鎖モンテカルロ法と階層ベイズモデルについて解説します。

事前分布から事後分布へ

4回中3回当たると、信じる度合いが0.7のほうへ動く

はじめに、ベイズ法の基本の式について解説します。

これまでの推定では、母数 $\theta$ を未知の定数と考えてきました。ベイズ法では、$\theta$ の不確かさを確率分布で表します。データを見る前の $\theta$ の分布を 事前分布（prior distribution）、データ $y$ を見たあとの分布を 事後分布（posterior distribution）と呼びます。

2つをつなぐのが、ベイズの定理です。

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

$\pi(\theta)$ が事前分布、$f(y \mid \theta)$ が尤度、$\pi(\theta \mid y)$ が事後分布です。分母はデータの周辺尤度で、$\theta$ を含まない定数です。そのため、次の一言に要約できます。

$$\text{事後分布} \propto \text{尤度} \times \text{事前分布}$$

天気予報の更新に例えられます。朝の予報が事前分布です。昼の雲の様子というデータを見て、午後の予報を修正します。これが事後分布です。そして、その事後分布は、次のデータを見るときの事前分布になります。

具体的に見てみましょう。くじの当たる確率 $\theta$ が 0.3, 0.5, 0.7 のどれかで、事前確率を 0.25, 0.5, 0.25 とします。4回引いて3回当たりました。

$\theta$	事前確率	尤度 $\theta^3(1 - \theta)$	事前 × 尤度	事後確率
0.3	0.25	0.0189	0.004725	0.077
0.5	0.50	0.0625	0.031250	0.506
0.7	0.25	0.1029	0.025725	0.417

表の見方はシンプルです。事前確率に尤度を掛けて、合計が1になるように割るだけです。0.3の確率が0.25から0.077に下がり、0.7の確率が0.25から0.417に上がりました。データが、信じる度合いを動かしたのです。

事後分布が得られたら、1つの値にまとめて報告します。

【事後分布のまとめ方】

• 事後平均： 二乗誤差の期待値を最小にする点推定
• 事後中央値： 絶対誤差の期待値を最小にする点推定
• 事後モード（MAP推定値）： 事後分布の密度が最大の点
• 信用区間： 事後確率が95%となる区間

信用区間（credible interval）は、「$\theta$ がこの区間に入る確率は95%」と、そのまま読める区間です。冒頭の信頼区間との違いは、ここにあります。

以上が、事前分布と事後分布についての解説です。

共役事前分布：手で計算できる組み合わせ

鍵と鍵穴の組み合わせなら、足し算でBeta(9,5)に更新できる

次に、事後分布が公式で求まる場合を説明します。

事後分布が事前分布と同じ分布族になるとき、その事前分布を 共役事前分布（conjugate prior）と呼びます。鍵と鍵穴がぴったり合う組み合わせのようなものです。パラメータを足し算で更新するだけで、事後分布が求まります。

代表的な組み合わせは次のとおりです。

データの分布	共役事前分布	事後分布
二項分布（$n$ 回中 $x$ 回）	$\text{Beta}(a, b)$	$\text{Beta}(a + x,\ b + n - x)$
ポアソン分布（$n$ 個の合計 $s$）	$\text{Ga}(a, b)$	$\text{Ga}(a + s,\ b + n)$
正規分布（分散 $\sigma^2$ は既知）	$N(m_0, \tau_0^2)$	正規分布

ガンマ分布は、$\theta^{a-1}e^{-b\theta}$ に比例する形（$b$ はレート）です。

1. ベータ分布と二項分布

なぜベータ分布になるのか、確かめましょう。尤度と事前分布の $\theta$ を含む部分を掛けます。

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

指数が足されるだけで、これは $\text{Beta}(a + x,\ b + n - x)$ の形です。$a$ と $b$ は「事前に見た当たりとはずれの回数」のように働きます。

例題：ベータ事前分布で当たる確率を推定する

事前分布を $\text{Beta}(2, 2)$ として、くじを10回引いたら7回当たりました。事後分布、事後平均、MAP推定値を求めましょう。

【解き方】

1. 事後分布： $\text{Beta}(2 + 7,\ 2 + 3) = \text{Beta}(9, 5)$
2. 事後平均： $\text{Beta}(a, b)$ の平均は $a/(a + b)$ なので、$9/14 \approx 0.643$
3. MAP推定値： モードは $(a - 1)/(a + b - 2)$ なので、$8/12 \approx 0.667$

答えは、事後分布 $\text{Beta}(9, 5)$、事後平均 約0.643、MAP推定値 約0.667です。ベータ分布の分位点を数値計算すると、95%信用区間は約 (0.386, 0.861)、$P(\theta > 0.5 \mid x) \approx 0.87$ です。

ここがポイントです。事後平均は、事前平均0.5と最尤推定値0.7の重みつき平均になっています。

$$\frac{9}{14} = \frac{4}{14} \times 0.5 + \frac{10}{14} \times 0.7$$

重みの4は事前分布の $a + b$、10はデータの数です。データが増えれば、事前分布の影響はだんだん薄れます。

2. 正規分布の平均

せっかくなので、別の例も見てみましょう。$y_1, \dots, y_n$ が $N(\mu, \sigma^2)$ に従い、$\sigma^2$ は既知で、事前分布を $\mu \sim N(m_0, \tau_0^2)$ とします。分散の逆数である 精度 を使うと、事後分布はきれいに書けます。

$$\mu \mid y \sim N\left(\frac{m_0/\tau_0^2 + n\bar{y}/\sigma^2}{1/\tau_0^2 + n/\sigma^2},\ \ \frac{1}{1/\tau_0^2 + n/\sigma^2}\right)$$

言葉で言うと、事後平均は事前平均とデータの平均を、それぞれの精度で重みづけた平均です。事後の精度は、2つの精度の和です。綱引きのように、確信の強いほうへ引っ張られます。

例えば、$m_0 = 60$、$\tau_0 = 5$、$\sigma = 10$、$n = 4$、$\bar{y} = 70$ とします。事前の精度は $1/25 = 0.04$、データの精度は $4/100 = 0.04$ で同じです。事後平均はちょうど真ん中の65、事後分散は $1/0.08 = 12.5$ です。

事前分布の選び方も確認しましょう。過去の知見があれば、それを反映した情報のある事前分布を使います。知見がなければ、一様分布 $\text{Beta}(1, 1)$ のように、なるべく影響の小さい 無情報事前分布 を使います。

試験のポイント：事後分布を求める問題では、$\theta$ を含む部分だけを掛けて、どの分布の形かを見抜きます。定数は最後まで無視してかまいません。

以上が、共役事前分布についての説明です。

MCMC：事後分布から乱数を引く

登りは必ず進み、下りは比の確率で進む

続いて、公式で解けないときの計算方法を解説します。

共役でない場合、ベイズの定理の分母の積分は、ふつう計算できません。母数が多いと、数値積分もお手上げです。ここで発想を変えます。事後分布の式が書けなくても、事後分布に従う乱数がたくさんあれば、その平均や分位点で事後分布を調べられます。

そこで登場するのが、マルコフ連鎖モンテカルロ法（MCMC法、Markov chain Monte Carlo）です。マルコフ連鎖の記事で、連鎖を長く走らせると定常分布に落ち着くことを学びました。それを逆に使います。定常分布が目標の事後分布になるようにマルコフ連鎖を設計して走らせ、その軌跡を乱数として使うのです。

その代表が MHアルゴリズム（メトロポリス・ヘイスティングス法、Metropolis-Hastings）です。目標の分布を $\pi(x)$、現在の値を $x$ とします。

【MHアルゴリズムの手順】

1. 提案分布 $q(y \mid x)$ から、候補 $y$ を引く
2. 採択確率 $\alpha(x, y)$ を計算する
3. 確率 $\alpha$ で $y$ に移動し、それ以外なら $x$ にとどまる
4. 1〜3を繰り返す

$$\alpha(x, y) = \min\left\{1,\ \frac{\pi(y)\, q(x \mid y)}{\pi(x)\, q(y \mid x)}\right\}$$

$q$ が左右対称なら、比は $\pi(y)/\pi(x)$ です。登山に例えるなら、目隠しで1歩踏み出す場所を決めて、登りなら必ず進み、下りならくじで進むかどうかを決めます。その結果、高いところほど長く滞在し、滞在時間が標高に比例するようになります。

大切なのは、$\pi$ が比でしか現れないことです。事後分布の定数である分母の積分は約分されて消えるので、「尤度 × 事前分布」だけで計算できます。

なぜ定常分布が $\pi$ になるのでしょうか？ $x$ から $y$ へ移る確率を $p(x, y) = q(y \mid x)\,\alpha(x, y)$ とすると、次の 詳細釣り合い条件 が成り立ちます。

$$\pi(x)\, p(x, y) = \min\{\pi(x)q(y \mid x),\ \pi(y)q(x \mid y)\} = \pi(y)\, p(y, x)$$

中央の式は $x$ と $y$ を入れ替えても同じです。両辺を $x$ について足すと、$\pi$ が定常分布であることが分かります。

具体的に見てみましょう。先ほどの $\text{Beta}(9, 5)$ を目標に、定数を除いた $x^8(1 - x)^4$ だけで動かします。現在が0.6で、対称な提案から候補が出たとします。

【1ステップの計算】

• 候補が0.7のとき： 比は約1.09なので、採択確率は1で必ず移動
• 候補が0.5のとき： 比は約0.57なので、確率0.57で移動
• 候補が0.9のとき： 比は約0.10なので、確率0.10で移動

標準偏差0.15の正規乱数で提案して20万回動かすと、平均は約0.642、標準偏差は約0.124となりました。真の値は0.643と0.124です。公式を使わなくても、事後分布が再現できています。

一方で、MCMCには実務上の注意があります。

【MCMCの注意点】

• 初期値の影響が残る最初の部分は捨てる（バーンイン）
• 連続する値には自己相関があるので、見た目の個数ほど情報がない
• 提案の幅が狭いと採択されても動きが遅く、広いと棄却が増える
• 初期値を変えた複数の連鎖が同じ分布に落ち着くかを確認する

試験のポイント：MHアルゴリズムの採択確率は比で決まるので、正規化定数が分からなくても使えます。棄却されたときは、「現在の値をもう1回記録する」点に注意します。

以上が、MCMCについての解説です。

ギブスサンプリング：1変数ずつ順番に更新する

飛車のように縦と横へ、1変数ずつ順番に動かす

さらに、MCMCのもう1つの代表を説明します。

母数が複数あると、全部を同時に提案するのは大変です。そこで、1変数ずつ順番に更新します。これが ギブスサンプリング（Gibbs sampling）です。

【ギブスサンプリングの手順（2変数の場合）】

1. 現在の $y^{(t)}$ を固定して、$x^{(t+1)}$ を $\pi(x \mid y^{(t)})$ から引く
2. 新しい $x^{(t+1)}$ を固定して、$y^{(t+1)}$ を $\pi(y \mid x^{(t+1)})$ から引く
3. 1と2を繰り返す

このように、他の変数すべてを固定したときの分布を 完全条件付き分布 と呼びます。縦に動いて、横に動いて、を繰り返す将棋の飛車のような動きです。

ギブスサンプリングは、MHアルゴリズムの特別な場合です。提案分布を完全条件付き分布にすると、採択確率の比は次のように1になります。

$$\frac{\pi(x', y)\, \pi(x \mid y)}{\pi(x, y)\, \pi(x' \mid y)} = \frac{\pi(x' \mid y)\, \pi(y)\, \pi(x \mid y)}{\pi(x \mid y)\, \pi(y)\, \pi(x' \mid y)} = 1$$

つまり、棄却が一度も起きません。その代わり、完全条件付き分布から乱数を引けることが条件です。

例えば、$(X, Y)$ が平均0、分散1、相関 $\rho$ の2変量正規分布なら、$X \mid Y = y \sim N(\rho y,\ 1 - \rho^2)$ です。$\rho = 0.8$ で今の $y$ が1なら、次の $x$ は $N(0.8,\ 0.36)$ から引きます。$x = 1.1$ が出たら、次の $y$ は $N(0.88,\ 0.36)$ から引きます。

せっかくなので、ベイズ推定の例も見てみましょう。正規分布の平均 $\mu$ と精度 $\phi = 1/\sigma^2$ がともに未知で、事前分布を $\mu \sim N(m_0, \tau_0^2)$、$\phi \sim \text{Ga}(a, b)$ と独立に置きます。同時の事後分布は名前のある形になりませんが、完全条件付き分布は名前のある形です。

$$\mu \mid \phi, y \sim N\left(\frac{m_0/\tau_0^2 + n\phi\bar{y}}{1/\tau_0^2 + n\phi},\ \frac{1}{1/\tau_0^2 + n\phi}\right), \qquad \phi \mid \mu, y \sim \text{Ga}\left(a + \frac{n}{2},\ b + \frac{1}{2}\sum_i (y_i - \mu)^2\right)$$

$\mu$ を見るときは $\phi$ を既知とみなし、共役事前分布の節の公式を使っています。$\phi$ のほうも同じ考えです。難しい同時分布を、簡単な1変数の問題に分けて解くのが、ギブスサンプリングの強みです。

ただし、変数どうしの相関が強いと、1回の動きが小さくなり、連鎖がなかなか全体を回りません。$\rho = 0.99$ なら条件つき分散は0.0199しかなく、細い谷を少しずつ進むことになります。

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

階層ベイズモデル：似たグループで情報を分け合う

口コミの少ない店ほど、全店の平均へ強く引き寄せられる

最後に、MCMCが特に力を発揮するモデルを見ていきましょう。

ラーメン店の口コミサイトを考えてみてください。A店は口コミ4件で平均75点、B店は口コミ100件で平均75点です。全店の平均は60点だとします。A店とB店の実力を、同じ75点と考えてよいでしょうか？ A店の4件は、たまたま好意的な客だったのかもしれません。

階層ベイズモデル（hierarchical Bayes model）は、事前分布の母数にもさらに事前分布を置き、グループの母数を2段階で表すモデルです。店 $j$ の $i$ 件目の点数を $y_{ij}$、店の実力を $\theta_j$ とすると、次のように書けます。

$$y_{ij} \mid \theta_j \sim N(\theta_j, \sigma^2), \qquad \theta_j \mid \mu, \tau^2 \sim N(\mu, \tau^2), \qquad (\mu, \tau^2) \sim \text{超事前分布}$$

1段目が店の中の点数のばらつき、2段目が店どうしの実力のばらつきです。3段目の、事前分布の母数の事前分布を 超事前分布（hyperprior）と呼びます。「どの店も、ラーメン店という同じ業界の出身」と考えるイメージです。

$\mu$、$\tau^2$、$\sigma^2$ を与えたとき、$\theta_j$ の事後平均は、共役事前分布の節の公式から次のようになります。

$$E[\theta_j \mid y] = w_j\, \bar{y}_j + (1 - w_j)\, \mu, \qquad w_j = \frac{n_j/\sigma^2}{n_j/\sigma^2 + 1/\tau^2}$$

$\sigma^2 = 100$、$\tau^2 = 25$、$\mu = 60$ として計算してみましょう。

店	件数 $n_j$	店の平均	重み $w_j$	事後平均
A店	4	75	0.50	67.5
C店	25	75	0.86	72.9
B店	100	75	0.96	74.4

件数の少ない店ほど、全店の平均に引き寄せられています。これを 縮小（shrinkage）と呼びます。店ごとに別々に推定すると少数のデータに振り回されます。全店を1つの平均で見ると店の差を無視します。階層モデルは、その中間で、データの量に応じて情報を分け合います。これを部分プーリングと呼びます。

実際には、$\mu$ と $\tau^2$ も未知です。$\theta_1, \dots, \theta_J, \mu, \tau^2$ と母数が多いので積分は手に負えませんが、完全条件付き分布は簡単な形になることが多く、ギブスサンプリングで計算できます。これが、MCMCと階層ベイズモデルの相性がよい理由です。

補足（ここは読み飛ばしてOK）： $\mu$ と $\tau^2$ に超事前分布を置かず、データから推定した値を代入して使う方法を、経験ベイズ法と呼びます。

試験のポイント：階層ベイズモデルでは、データの少ないグループの推定値ほど、全体の平均へ強く縮小されます。

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

細かい式を忘れてしまっても問題ありません。ここでは、「事後分布は尤度 × 事前分布」「解けないときはMCMCで乱数を引く」「階層モデルはグループで情報を分け合う」という点だけ押さえておきましょう。

今回のまとめ

尤度×事前で更新し、解けなければMCMCで乱数を引く

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

今回は、母数の不確かさを確率分布で表し、データで更新するベイズ法について紹介しました。やっていることは、天気予報の更新と同じです。計算の難しさはMCMCが引き受けてくれるので、複雑なモデルでも「事前分布と尤度を書く」ことに集中できます。

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

振り返り

• 事後分布は「尤度 × 事前分布」に比例する。信用区間は「母数が入る確率」とそのまま読める。
• 共役事前分布なら、パラメータの足し算で事後分布が求まる。事後平均は事前平均とデータの重みつき平均。
• MCMCは、定常分布が事後分布になるマルコフ連鎖を走らせて、乱数を得る方法。
• MHアルゴリズムは比で採択を決めるので、正規化定数がいらない。登りは必ず進み、下りはくじで決める。
• ギブスサンプリングは完全条件付き分布から順に引く方法で、棄却がない。
• 階層ベイズモデルは、データの少ないグループを全体の平均へ縮小する。

キーワードのおさらい

• 事前分布: データを見る前の母数の分布
• 事後分布: データを見たあとの母数の分布。尤度 × 事前分布に比例する
• 共役事前分布: 事後分布が事前分布と同じ分布族になる事前分布
• マルコフ連鎖モンテカルロ（MCMC）法: 目標の分布を定常分布にもつマルコフ連鎖で乱数を得る方法
• MH（Metropolis-Hastings）アルゴリズム: 提案された候補を採択確率で受け入れるMCMCの方法
• ギブスサンプリング: 完全条件付き分布から1変数ずつ順に乱数を引くMCMCの方法
• 階層ベイズモデル: 事前分布の母数にも事前分布を置く、多段階のモデル

さて、MCMCでは、乱数がたくさん引ければ分布が分かることを見ました。でも、そもそもコンピュータはどうやって乱数を作っているのでしょうか？ 乱数で計算する方法は、ほかにもあるのでしょうか？

次回は、乱数の生成法と「モンテカルロ法」「ブートストラップ法」などの計算手法について紹介します。

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