密度が計算できてもサンプリングできない：MCMCを自作した記録【番外編】

はじめに

第33回でベイズ法を扱ったとき、MCMC（Markov Chain Monte Carlo＝マルコフ連鎖モンテカルロ法）の考え方を追いました。ただ、あのときはアルゴリズムを読んだだけでした。

そこで練習用のプロジェクトを作り、メトロポリス法の中核だけを自分で書くという課題をやりました。目標分布・診断図・簡易チェックは用意済みで、自分が書くのは受理判定のループだけ、という形です。

書き終えて分かったのは、自分がいちばん誤解していたのはアルゴリズムではなく「MCMCが何のための道具か」だったということです。とくに次の2つを混同していました。

• 確率密度を点ごとに計算できること
• その分布から標本を生成できること

この2つはまったく別です。そして、この区別が付いていなかったので「解析的に解ける正規分布に、なぜわざわざMCMCを使うのか」がずっと引っかかっていました。

なお用語について1つ。第33回では「受容」「受容率」と書きましたが、この記事では実装のコードに合わせて「受理」「受理率」で統一します。同じものです。

この記事は、その引っかかりを実装と数値で解いた記録です。ベイズの理屈やギブスサンプリング、$\hat{R}$ による収束判定は第33回、逆関数法・棄却法・重点サンプリングそのものは第34回で扱ったので、ここでは重複を避けて「使い分け」に絞ります。

結論を先に3つ

① MCMCは確率密度関数を発見する方法ではありません。 目標分布に従う、互いに相関のある標本を生成する方法です。密度の式は最初から手元にあります（なければMCMCも動きません）。

② 目標分布から独立な標本を直接生成できるなら、通常は直接モンテカルロのほうが簡単で効率的です。 後で数値を出しますが、同じ標本数でも実質の情報量が20分の1近く違いました。

③ MCMCの利点は、正規化定数が分からなくても、正規化前の密度の相対値を評価できれば使えることです。 ベイズの事後分布がまさにその形をしています。（ただしこの性質自体はMCMC固有ではありません。後で触れます。）

以下、この3点を実装と数値で確かめていきます。

この記事で扱う用語

密度を評価できることと、サンプリングできることは別

まず、この記事の出発点になった混同をはっきりさせます。

「密度を評価できる」とは、$x$ を1つ渡せば $p(x)$（または正規化前の値）を計算して返せる、ということです。 関数に代入するだけなので、式が書けていれば必ずできます。

「サンプリングできる」とは、その分布に従う値を1つ生成できる、ということです。 これは代入では作れません。乱数の生成器が必要です。

ベイズの事後分布は、この表の右下に落ちます。 尤度と事前分布を掛ければ任意のパラメータ値での「もっともらしさ」は計算できます。しかし「事後分布から1個引く」ための専用の乱数生成器は、一般には存在しません。

MCMCはこの穴を埋める道具です。評価はできるがサンプリングできない分布から、標本を作る。

通常のモンテカルロ法・重点サンプリング・MCMCの違い

3つを並べて整理します。重点サンプリングはMCMCの一種ではありません。 別系統の手法です。

$Z$ の行に注意が必要です。「正規化定数が不要」はMCMC固有の性質ではありません。 自己正規化した重点サンプリング（重みを合計1に直す形）も $Z$ を必要としません。だから両者の本当の分かれ目は $Z$ ではなく、重みが壊れるかどうかでした。

同じ推定対象を3手法で

目標分布を $N(2, 1.5^2)$ とし、$E[X]$（真値 2.0）を推定します。標本数はすべて 9000 に揃えました。

②の標準誤差は自己正規化のため近似値です。中身は「重み付き分散 ÷ 重みのESS」の平方根で、デルタ法の簡略形にあたります。①③は通常の標準誤差です。

この1次元の問題では ① が圧勝です。 同じ9000標本でも、実質の情報量が 9000 対 4571 対 526。目標分布から直接引けるなら、MCMCを使う理由はありません。

では差はどこで出るのか

同じ形の問題を $d$ 次元にします。目標を $N(2\mathbf{1}, 1.5^2 I)$、重点サンプリングの提案を $N(0, 3^2 I)$ とすると、重みは各座標の比の積になります。

（重点サンプリング側の分母は2万、MCMC側は burn-in 2000 を落とした1万8000です。）

$d = 50$ では、2万個引いた重みのうち1個が全体の73%を占めます。 残りの 19,999 個を全部合わせても 27% ほどで、重みのESSは実質2個分（$0.0001 \times 20000$）しかありません。重点サンプリングは次元とともにこう壊れます。一方MCMCのESS比はゆるやかに下がるだけです。

（MCMCの側は提案幅を $2.4/\sqrt{d}$ で調整しています。表の $d=1$ の値 0.1770（17.70%）が、さきほどの3手法比較で出た 5.8% と違うのは提案幅が 2.4 と 1.0 で違うためです。
なお理論的な目安は $2.38\,\sigma/\sqrt{d}$ で、$\sigma = 1.5$ なら $d=1$ で約 3.6 になります。ここでは $\sigma$ を掛けずに使ったので最適よりやや狭く、17.70% は後で出てくる最良の 22.02% に届いていません。なお2つは鎖の長さが違うので、厳密な比較ではなく目安の比較です。）

これが使い分けの分かれ目でした。 重点サンプリングとMCMCの違いは「$Z$ が要るかどうか」ではなく、「重みで補正するか、受理判定で補正するか」の違いで、後者のほうが次元に強い。

なぜ既知の正規分布を教材にするのか

練習用プロジェクトの目標分布はこうなっています。

Target: Normal(mean=2.0, std=1.5)
log_target(x) = -0.5 * ((x - 2.0) / 1.5)^2

正規分布なら直接サンプリングできますし、平均も標準偏差も解析的に分かっています。 つまり実用上、この分布にMCMCを使う必要はまったくありません。

それでも教材にする理由は1つです。答えが分かっているので、実装が正しいかを確かめられる。

実際に、用意されていた簡易チェック（check_sampler.py）は次を検査します。

標本平均:   1.9811（理論値 2.0000）
標本標準偏差: 1.5147（理論値 1.5000）
受理率:     79.9%
OK: 平均
OK: 標準偏差
OK: 受理率
OK: チェーンが動いている

すべての簡易チェックを通過しました。

答えの無い分布で書いたら、この4行が出せません。「正しく動いているか分からないアルゴリズム」を書いてしまうのがいちばん怖いので、最初は答えのある分布で書くというのが教材の意図でした。

$p(x)$、$S(x)$、$\log p(x)$ の対応

ここで用語を1つ整理しておきます。物理由来の書き方で、目標分布を次の形に書くことがあります。

この $S(x)$ をエネルギー関数と呼びます。両辺の対数を取ると、

つまり $S(x) = -\log p(x) + \text{定数}$ と置いただけで、特別な仮定は入っていません。正規分布なら

です。$Z$ が正規化定数で、$\int p = 1$ にするための割り算の分母です。

3枚とも同じ分布を描いています。左は縦の縮尺が $Z$ 倍違うだけ、右は上下に $\log Z$ だけ平行移動しているだけです。

数値で確かめます。

密度では「比」、エネルギーと対数密度では「差」

ここが実装で間違えやすいところです。受理判定に使うのは比です。

$Z$ が約分されて消えます。 そして指数の中は $S$ の差です。3通りの移動で確かめました。

読み取れることが3つあります。

• 密度が上がる移動（$S$ が下がる移動）は比が1以上になり、必ず受理されます。 1行目と3行目です
• 密度が下がる移動でも、比の分だけの確率で受理されます。2行目は 22.3% で受理
• 右端の「差」の列は受理判定に使いません。 $p(x_{\text{new}}) - p(x)$ は $Z$ が約分されないので、正規化定数が必要になってしまいます

確率密度そのものでは「比」を取り、エネルギーまたは対数密度では「差」を取る。 これが対応関係です。実装では対数で持って差を取るので、桁溢れも起きません。

ランダムウォーク・メトロポリス法：自分で書いた部分

私が書いたのは次の5ステップだけです。目標分布・診断図・チェックは用意されていました。

current = initial_state
current_log_density = log_target(current)

for i in range(n_samples):
    dx = rng.normal(loc=0.0, scale=proposal_std)   # 対称な正規ランダムウォーク
    candidate = current + dx
    candidate_log_density = log_target(candidate)

    log_alpha = candidate_log_density - current_log_density   # 対数受理比

    r = rng.random()
    if min(0.0, log_alpha) > np.log(r):             # log(U) < min(0, log_alpha)
        accepted_count += 1
        current = candidate
        current_log_density = candidate_log_density

    samples[i] = current      # 受理・棄却にかかわらず、現在の状態を保存する

書いてみて引っかかった点が3つありました。

① 棄却したときも標本に記録する。 ここを「受理したときだけ記録」にすると別の分布になります。棄却は「同じ場所にもう1歩ぶんの滞在時間を与える」という意味を持っているので、記録しないと滞在時間の情報が消えます。

② 比較は対数で行う。 $\alpha = \exp(\log\alpha)$ を計算して $U < \alpha$ と比べても数学的には同じですが、$\log\alpha$ が大きな負の数のとき $\exp$ が 0 に潰れます。$\log U < \min(0, \log\alpha)$ の形なら潰れません。

③ min(0, ...) は要るのか。 受理確率は $\min(1, \alpha)$ なので、対数では $\min(0, \log\alpha)$ です。$\log U$ は必ず負なので、$\log\alpha > 0$ のときは $\min$ を取らなくても常に受理されます。つまり結果は同じですが、式が受理確率の定義と1対1に対応しているので付けたままにしました。

対称な提案分布だから Hastings 補正が消えている

一般のメトロポリス–ヘイスティングス法の受理確率は次の形です。

右側の分数が Hastings 補正です。今回の提案分布は「現在位置を中心とする正規分布」なので対称で、

が成り立ちます。だから補正が1になって消え、古典的なメトロポリス法（Metropolis et al. 1953）の形になります。

非対称な提案を使うときは補正が必要です。 たとえば正の値しか取らないパラメータに対して対数正規の提案を使う場合や、現在位置に応じて提案の幅を変える場合です。補正を忘れると別の分布に収束します。しかも「動いてはいる」ので、気づきにくい。

なお、対称性は必須条件ではありません。目標分布を不変分布にするための十分条件は詳細釣り合いで、「対称な提案＋メトロポリスの受理則」はそれを満たす組み合わせの一例にすぎません。 詳細釣り合いを満たすのは提案分布単体ではなく、提案と受理則を合わせた遷移全体です。詳細釣り合いそのものは第33回で扱ったので、ここでは深追いしません。

提案分布を変えても同じ分布に収束するのはなぜか

これも引っかかった点でした。提案分布はまったく違うものを使えるのに、なぜ同じ目標分布に行き着くのか。

理由は、提案分布は「どこを次に見るか」しか決めておらず、「そこに留まる割合」は受理判定が決めているからです。受理確率に $p$ の比が入っているので、最終的な滞在時間の比が $p$ の比に一致します。提案分布は探索の効率だけを変えます。

ただし無条件ではありません。次の条件が必要です。

ランダムウォーク・メトロポリス法では、提案が全域に正の密度を持つ正規分布なので既約性と非周期性は自動的に満たされます。だから実務で気をつけるのは1行目、受理確率を正しく書くことだけになります。

診断図の読み方

固定シードで実行した結果を見ます。設定は --samples 10000 --burn-in 1000 --proposal-std 1.0 --seed 42（初期値 0）です。

4枚それぞれの読み方を書きます。

① ヒストグラムと理論密度

見るのは「形が合っているか」だけです。 ここが合っていなければ実装が間違っています。合っていても収束の証明にはなりません（後述）。

② トレースプロット

縦軸が状態、横軸が反復回数です。 見るポイントは3つ。

• 初期値からの立ち上がり：0 から始めて、最初の20歩は 0.305, 1.055, 1.055, 1.183, 1.166, 2.046, … と目標付近へ上がっていきます。同じ値が2回続いているところ（1.055 が2回）は棄却された歩です
• 同じ場所に長く留まっていないか：水平な線が長く続くのは棄却が連続しているサインです
• 帯の太さが一定か：途中で急に幅が変わるなら、分布の別の領域に移った可能性があります

赤い縦破線が burn-in の終わりですが、これはグラフから自動検出した値ではありません。 コマンドラインで --burn-in 1000 と指定した値をそのまま描いているだけです。この実装に自動判定は入っていません。

③ 自己相関

MCMC標本は独立ではありません。 ここが通常のモンテカルロ法との最大の実務的な違いです。

lag 1 の自己相関は 0.8779。1歩前の自分と 88% 相関しています。0.05 を下回るのは lag 26 で、26歩離れてようやく「ほぼ無関係」になります。

これを1つの数にまとめたのがESSです。

$\rho_k$ が lag $k$ の自己相関です。自己相関が全部 0 なら ESS $= n$。相関が強いほど小さくなります。今回は ESS = 526、つまり9000標本が実質526個分の情報しかありませんでした。

（無限和は打ち切る必要があります。この記事では Geyer (1992) の「連続する2つの自己相関の和が負になったら止める」という方式を使っています。）

ESSが分かると平均のズレが解釈できます。 モンテカルロ標準誤差は標準偏差 ÷ $\sqrt{\text{ESS}}$ なので、$1.5837 / \sqrt{526} = 0.0691$。実測の平均 2.0490 は理論値から 0.0490 離れていますが、これは 0.71 標準誤差ぶんです。偏りではなく誤差の範囲でした。

④ ランニング平均

先頭から $k$ 個までの平均を $k$ の関数として描いたものです。理論平均 2.0 に近づいていく様子が見えます。

ただし、ここが最も誤解しやすい図です。 ランニング平均が横ばいになっても、それは収束の証明になりません。実測すると、$\pm 0.05$ の帯に

• はじめて入るのは 144 標本目
• それ以降ずっと帯の中に留まるのは 8451 標本目

でした。144 標本目のグラフを見て「安定した」と判断すると間違えます。横ばいに見えるのは、単に分母が大きくなって新しい標本の影響が薄まっているからでもあります。

受理率

今回は 79.5% でした。4枚のパネルには出ませんが、図全体の見出しに表示されています。この数字の読み方は次の節でまとめます。

実務での収束判定

この記事の診断（1本の鎖のトレース・自己相関・ランニング平均）は、どれも「明らかに壊れている」ことを見つける道具です。 「収束した」ことの証明にはなりません。

実務では次を併用します。

$\hat{R}$ は第33回で扱いました。今回のプロジェクトは1本の鎖しか走らせないので、原理的に $\hat{R}$ が計算できません。 そこは教材の範囲外だと理解しておく必要があります。

提案幅を変えると何が変わるのか

同じ目標分布・同じシード・同じ標本数で、提案分布の標準偏差だけを変えます。

（上の図の見出しには「ESSが最大なのは 3.0」と書いてありますが、これは1万標本での比較です。あとでもっと細かく探すと 3.5 が最良になります。）

受理率が最も高い 0.2 が、いちばん効率が悪いです。 ESSが 32 しかありません。9000歩使って実質32個分。

理由はトレースを見れば分かります。提案幅 0.2 では、ほとんどの候補が受理されるかわりに、1歩で 0.2 程度しか動きません。 ゆっくり漂うだけで、2000歩でも分布の全域を回りきれない。「受理された」のと「進んだ」のは別です。

逆に提案幅 10.0 では、82% の反復が棄却されます。棄却は「同じ場所で足踏み」なので、これも効率を落とします。

つまり効率には山があります。 もっと細かく探しました（10万標本）。

この1次元の問題では、提案幅 3.5（受理率 45.0%）が最良でした。 目標分布の標準偏差 1.5 に対して 2.33 倍です。

高次元では漸近的に最適な受理率が 0.234 であることが知られています（Roberts, Gelman & Gilks 1997）。この実測はそれよりずっと高い受理率で最良になりました。 ただしこれは理論と食い違っているわけではありません。1次元での最適な受理率は約 0.44 とされており（Roberts & Rosenthal 2001）、実測の 45.0% はそちらとよく一致しています。

つまり「受理率を 0.234 に合わせる」は高次元での目安で、低次元ではもっと高いほうが良い。 どちらにしても受理率だけを見て調整せず、ESSで比べるのが確実です。

極限の分布は変わらない

「提案幅で結果が変わるなら、正しい分布に収束していないのでは」と思うかもしれません。標本数を10倍にして確かめました。

どの提案幅でも平均が 2.0 に、標準偏差が 1.5 に近づきます。 提案分布は極限の目標分布を変えません。変えるのは有限標本での探索効率だけです。

さらに、1万標本のときの平均のズレも誤差として説明できます。

どれも1標準誤差以内です。 提案幅 0.2 の 2.2024 は偏りではなく、ESSが32しかないための誤差でした。平均のズレはESSを見ないと解釈できません。

密度は評価できるが直接サンプリングできない例

ここまでは正規分布でしたが、そこにMCMCを使う実用上の理由はありませんでした。実際に必要になる形を1つ見ます。ロジスティック回帰のベイズ事後分布です。

設定

説明変数 $x$ と二値の応答 $y$ を20組用意しました。

x = -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.0 0.5 1.0 1.5 2.0 2.5 -1.2 0.7 1.8 -0.3 1.1 -2.2 0.2 2.2
y =  0    0    0    0    1   0   1   1   1   1   1   1    0   0   1    1   1    0   0   1

モデルは次の形です。$\theta = (\text{intercept}, \text{slope})$ の2次元です。

事前分布は intercept・slope ともに $N(0, 2.5^2)$ の独立としました。これは固有な（proper な）事前分布です。

事前分布を固有にしたのには理由があります。ロジスティック回帰では完全分離（ある閾値で $y$ が完全に 0 と 1 に分かれる状態）が起きると、平坦な事前分布のもとで事後分布が発散して不適切になります。今回のデータは $x$ の順に並べると

0, 0, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1

で、どの閾値でも完全には分かれません。そのうえ固有な事前分布を置いているので、事後分布は二重に安全に proper です。

密度は必ず計算できる

対数事後密度は「対数尤度 + 対数事前」です。$\theta$ を渡せば必ず値が返ります。

しかし「この事後分布から $\theta$ を1つ引く」ための専用の乱数生成器は存在しません。 正規分布やベータ分布のようには作れない。ここが表の右下でした。

正規化定数は二重積分になる

正規化定数は

です。2次元なら格子で数値積分できます。ただしここで1つ落とし穴を踏みました。

実装の log_prior は事前分布の正規化定数を省いています（受理判定では消えるので不要だからです）。1成分あたり $-\frac{1}{2}\log(2\pi \cdot 2.5^2) = -1.835229$、intercept と slope の2成分で $-\log(2\pi \cdot 2.5^2) = -3.670459$ です。そのまま格子で積分すると $\log = -6.341098$、値 $= 1.762367 \times 10^{-3}$ になりますが、これは上の式の $Z$ ではありません。この $-3.670459$ を足すと

で、これが上の式が意味する値（ベイズの周辺尤度）です。2つは $2\pi \cdot 2.5^2 = 39.2699$ 倍違います。
なお厳密には、これも格子の範囲（intercept は $[-6,6]$、slope は $[-6,10]$）で打ち切った数値積分の値です。

受理判定では消える定数が、正規化定数そのものを求めるときには消えない。 この記事の主題（比では $Z$ が消える）の裏返しで、自分で1回踏みました。

ただしMCMCはこの値を一切使いません。 そしてパラメータが10個あれば10重積分になり、この方法は使えなくなります（第34回で見たグリッド法の次元の壁と同じ話です）。

2次元のランダムウォーク・メトロポリス法

1次元のときと同じ受理判定を、2次元にしただけです。提案は各成分に独立な正規分布（標準偏差 0.55）、6万反復、burn-in 1万、初期値 $(0, 0)$。

受理率は 62.9% でした。

格子による数値積分と照合しました。

一致しています。 2次元だから格子と比べられるので、これも「答えのある問題で実装を確かめる」の一例です。

標本があると何ができるか

事後分布の標本が手に入ると、知りたい量を標本に対して計算して数えるだけで済みます。

• 傾きが正である事後確率：$P(\text{slope} > 0) = 0.9999$。標本のうち slope が正だったものの割合を数えただけ
• $x = 1.0$ での予測確率：各標本の $\theta$ をシグモイド関数に通して平均すると 0.8655、95%区間は [0.6095, 0.9889]

右のパネルの水色の束がこれです。係数の標本1個ごとに予測曲線が1本引けるので、予測の不確実性がそのまま束の太さになります。 点推定1つでは出てこない情報です。

そしてこの計算のどこにも正規化定数 $Z$ は出てきません。 受理判定が比しか使わないので、$Z$ を知らないまま標本が得られ、標本から期待値も確率も出せます。

ただし $Z$ そのものは出ません。 通常のMCMCは正規化定数を推定しません。$Z$（ベイズでは周辺尤度）が欲しい場合は、ブリッジサンプリングや熱力学的積分といった別の工夫が必要です。

通常のモンテカルロ法で十分な例

逆の方向も見ておきます。球の体積を求める問題です。

球の条件は

で、立方体 $[-r, r]^3$ から一様に点を打って、球の中に入った割合を数えます。

3次元ではMCMCを使う必要がありません。 理由を正確に書きます。

「中心に確率が集中していないから」ではありません。 正しい理由は次の2つです。

1. 基準となる一様分布から独立標本を直接生成できる（rng.uniform で済む）
2. 球内に入る事象が希少ではない（命中率 52.4%）

このどちらかが崩れると話が変わります。

高次元では崩れる

球と超立方体の体積比は次元とともに急減します（第34回でも同じ比率を見ました）。

$d = 20$ では 4000万回引いて1個です。上の条件2（希少でない）が崩れています。

ただし「単純なMCMCに替えれば自動的に体積が求まる」わけではありません。 ここは私が最初に誤解しかけた点です。通常のMCMCは目標分布の標本を作るだけで、正規化定数を出しません。そして球の体積はまさに正規化定数の側にある量です。高次元の体積を求めるには、重点サンプリング・逐次的に領域を分割する方法・熱力学的積分といった、別の枠組みが必要になります。

使い分けの判断フロー

ここまでを1本の流れにまとめます。

目標分布から独立標本を直接生成できる？
  Yes → 通常のモンテカルロ法
  No  → 良い提案分布と安定した重みを作れる？
           Yes → 重点サンプリング／棄却サンプリングを検討
           No  → 正規化前の目標密度を評価できる？
                    Yes → MCMC を検討
                    No  → 通常の MCMC も困難

別枠の注意が1つあります。 上のフローの1行目が Yes でも、通常のモンテカルロ法が実用にならない場合があります。知りたい事象がまれなときです。

たとえば「確率 $10^{-8}$ の事象の確率を推定したい」なら、直接生成できても素朴に数えるだけでは1億回引いて1回しか当たりません。この場合は重点サンプリングなどで「まれな事象を起きやすくして重みで割り戻す」必要があります。第34回で $10^{-23}$ の確率を10万回で推定した例がこれです。

つまり判断は2軸で見るのが正確です。

技術的に混同しやすい点のまとめ

自分が実際に混同していた点、あるいは書きながら怖くなって確認した点を並べます。

参考文献

一次資料と標準的な教科書を挙げます。

再現方法

この記事の数値と図は、すべて固定シードで再現できます。スクリプトは site/scripts/mcmc-column/ に置いてあります。

mcmc_common.py は練習用プロジェクトのパスを絶対パスで直書きしているので、別の環境で動かすときはそこを書き換えてください。

主な設定は次のとおりです。

• 目標分布：$N(2.0, 1.5^2)$。target.py の log_target が返すのは $-\frac{1}{2}\left(\frac{x-2.0}{1.5}\right)^2$（正規化定数を省いた対数密度）
• 診断図：--samples 10000 --burn-in 1000 --proposal-std 1.0 --seed 42（初期値 0）
• 簡易チェック：n_samples=30000, burn_in=3000, proposal_std=1.0, initial_state=-5.0, seed=2026
• ロジスティック回帰：提案の標準偏差 0.55、6万反復、burn-in 1万、seed 20260816、事前分布 $N(0, 2.5^2)$

この記事の位置づけ

連載の中での位置を書いておきます。

準1級の試験対策としては、この記事の内容は「用語と使いどころ」のレベルで十分です。 MCMCを実装させる問題は出ません。ただ、自分で書いてみると「棄却したときも記録する」「対数で比較する」「受理率が高いほど良いのではない」といった、読むだけでは通り過ぎる点が体に入りました。

この連載の全体像は統計検定準1級・独学連載のまとめにあります。