確率・分布●●●●●

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

Markov Chain Monte Carlo ・ まるこふれんさもんてかるろほう

マルコフ連鎖で事後分布からサンプルを生成する計算手法。ベイズ実用化の要

✦ 30秒で試す — 航海歩幅が小さすぎても大きすぎても、山の間を渡れない歩幅が小さすぎて片方の山から出られないとき、受理率は?予想してみる

概要

マルコフ連鎖モンテカルロ法(MCMC)は、数式のままでは扱えない複雑な確率分布から、サンプル(乱数)を生成するための計算手法です。名前のとおり2つの部品の合体で、「狙った分布が定常分布になるように設計したマルコフ連鎖」を走らせ、そこから得たサンプルの集まりで分布を近似するモンテカルロ法を実行します。連鎖を長く走らせるほど、訪れた点の集まりが狙った分布の形を写し取っていきます。

MCMCが最も活躍する舞台はベイズ推論です。ベイズ推論の答えである事後分布は、理屈のうえでは事前分布と尤度の掛け算で決まりますが、確率として正規化するための分母(あらゆるパラメータ値にわたる積分)が、現実的なモデルではまず計算できません。MCMCの妙は、この計算できない分母を知らないまま、事後分布からサンプルを引き出せることです。

サンプルさえ手に入れば、あとは集計するだけです。平均を取れば事後平均、分位数を取れば信用区間、ヒストグラムを描けば分布の形。「積分を解く」という数学の問題が、「サンプルを集計する」というプログラミングの問題に置き換わる — これがMCMCの提供する変換です。

なぜ生まれたか

MCMCの原型は統計学ではなく物理学から生まれました。1953年、ロスアラモス研究所のメトロポリスらは、多数の粒子からなる系の平衡状態を計算する問題に直面していました。取りうる配置の数が天文学的で、まともに全配置を足し合わせることは不可能です。彼らは「配置を少しランダムに動かし、エネルギー的に有利なら受理、不利でも一定の確率で受理する」という単純な反復で、物理的に正しい分布どおりに配置をサンプリングできることを示しました。これがメトロポリス法で、世界初期の電子計算機MANIACの上で動きました。

一方、ベイズ統計は長らく「理論は美しいが計算できない」学問でした。事後分布の正規化定数が積分できるのは、共役事前分布が使えるような教科書的に単純なモデルに限られ、現実的な複雑さのモデルは手が出せなかったのです。1970年にヘイスティングスがメトロポリス法を統計の文脈に一般化し、画像処理の分野で1984年にジェマン兄弟が提案したギブスサンプリングを、1990年にゲルファンドとスミスがベイズ計算の一般的な道具として示したことで、状況は一変します。「どんなモデルでも、事後分布に比例する量さえ計算できればサンプリングできる」— この汎用性が、頻度論とベイズの長い論争の実務的な足かせだった計算問題を取り除き、1990年代以降のベイズ統計の復権を直接引き起こしました。

詳細

メトロポリス・ヘイスティングス法 — 受理と棄却の繰り返し

MCMCの最も基本的なアルゴリズムがメトロポリス・ヘイスティングス法です。現在の点の近くに候補の点をランダムに提案し、事後分布の密度(正確にはそれに比例する値)を現在の点と比べます。候補のほうが密度が高ければ移動し、低ければ「密度の比」を確率として移動するかどうかを乱数で決めます。この一歩を何万回も繰り返すだけです。

はい

いいえ

今の点から出発

提案分布から
近くの候補点を引く

比 r を計算
分子: 候補の密度 × 戻る提案の確率
分母: 今の密度 × 行く提案の確率

一様乱数 u が
r 未満か

受理して
候補点へ移動

棄却して
今の点に留まる

今いる点を
1標本として記録

正確に書くと、受理の判定に使う比 r は「候補の密度 × 候補から今の点へ戻る提案が出る確率」を「今の密度 × 今の点から候補への提案が出る確率」で割ったもので、受理確率は r と1の小さいほうです。提案確率の比(ヘイスティングスの補正)は、提案に偏りがあるときにそれを打ち消す役目をします。たとえば「右へ提案しやすい」提案分布を使うなら、右への移動を少し割り引かないと連鎖が右へ偏ってしまうからです。現在の点を中心にした正規分布のように、行きと戻りの提案確率が等しい左右対称な提案なら、この補正は1になって消え、密度の比だけで判定できます。これがヘイスティングスの一般化以前の、もともとのメトロポリス法です。

鍵は、判定に使うのが密度の比であることです。比を取ると、計算不能だった正規化定数(分母)が分子・分母で打ち消し合って消えます。だから「事後分布に比例する量 = 事前分布 × 尤度」さえ計算できれば連鎖を回せる。密度の高い領域には長く滞在し、低い領域はたまにしか訪れないため、記録された点の集まりが事後分布の形を再現します。マルコフ連鎖の言葉で言えば、この受理・棄却の規則は「どの2点の間でも、行きの流れと戻りの流れが釣り合う」(詳細釣り合い)ように遷移確率を設計したことに相当し、その結果として目標の分布が定常分布になります。棄却されたときに「同じ点をもう1回記録する」のを省くと、この釣り合いが崩れて分布がずれる点にも注意が必要です。

下の部品で、2次元の分布の上を歩き手に歩かせてみてください。歩幅(提案分布の広さ)を変えると、受理率と「分布全体への広がりやすさ」がトレードオフになっていることが見えてきます。小さすぎる歩幅は受理率こそ高いものの一向に進まず、大きすぎる歩幅は遠くへ飛びすぎてほとんど棄却されます。

⚡ 体験: メトロポリス法で分布の上を歩く
目標
歩幅
目標分布と歩き手0 歩トレース(x 座標)burn-in-4-2024歩数 →標本
歩数0受理率—x>0 にいた割合—目標45%

金の歩き手は分布のすそ野の外れから出発します(つかんで置き直せます)。「1歩」を押すと、近くに候補を出し、密度の比 r で受理・棄却を決めます。

等高線 = 目標分布(正規化定数は知らないまま)。金 = 受理された足跡、朱 = 棄却された候補、白い点の雲 = 集めた標本

主要なバリエーション

ギブスサンプリングは、パラメータが複数あるとき「他をすべて固定した1つのパラメータの条件付き分布」から順番にサンプリングする方式です。各条件付き分布が既知の形になるモデル(階層ベイズモデルの多くがそうです)では提案の棄却がなく効率的で、初期のベイズソフトBUGSの中核でした。ただしパラメータ同士の相関が強いと、1軸ずつしか動けないため細い尾根に沿ってジグザグにしか進めず、効率が落ちます。

**ハミルトニアン・モンテカルロ法(HMC)**は、事後密度の勾配(傾き)情報を使い、物理の運動方程式を模してパラメータ空間を滑らかに長距離移動する方式です。ランダムに一歩ずつよろよろ歩くメトロポリス法に比べ、パラメータ数が多い高次元の事後分布でも効率よく探索できます。その自動調整版であるNUTSが、Stan や PyMC といった確率的プログラミング言語の標準エンジンです。現代の実務では、モデルを記述すればサンプラーの選択と調整はツールに任せる、という分業が一般的になっています。

収束の診断 — 信じてよいサンプルか

MCMCには理論保証と実務のギャップがあります。定常分布への収束は「無限に走らせれば」の保証であり、有限のサンプルが事後分布を十分代表しているかは別途確かめる必要があります。実務の定石は、出発点をばらばらにした複数の連鎖(チェーン)を走らせ、それらが同じ分布に混ざり合ったかを見ることです。

よく混ざった連鎖 — 歩幅が適切(R-hat = 1.00)-404混ざっていない連鎖 — 歩幅が小さすぎる(R-hat = 2.36)-404反復回数 →burn-in標本 vs 目標反復回数 →山の比率が狂うθθ
トレースプロットで収束を診断する — 出発点の違う3本の連鎖が同じ帯に混ざれば合格。灰色の帯は捨てる burn-in、右は burn-in 後の標本のヒストグラムと目標の分布(破線)

まず、連鎖の序盤は出発点の影響が残っているため、バーンイン(ウォームアップ)として捨てます。次に、複数の連鎖の「連鎖どうしのばらつき」と「連鎖の中のばらつき」を比べるR-hat統計量で、混ざり具合を数値で確認します(1.01以下が目安)。図の下段のように、各連鎖が別々の山に閉じ込められていると、1本1本は落ち着いて見えても R-hat が1を大きく超えます。また、マルコフ連鎖のサンプルは隣どうしが相関しているため、独立なサンプル何個分の情報量に相当するかを有効サンプルサイズで見積もります。トレースプロット(サンプル値の時系列グラフ)が特定の領域に長く停滞していないかを目で確かめることも基本です。

注意したいのは、診断は「うまくいっていない証拠」を見つける道具であって、「うまくいっている証明」ではないことです。すべての連鎖が同じ1つの山に閉じ込められ、別の山の存在に誰も気づかなければ、R-hat は1のまま正常に見えます。診断を省いたり過信したりして収束していない連鎖を集計すると、もっともらしい数字の形をした誤答が得られます。MCMCの最大の落とし穴は、失敗しても一見それらしい結果が出てしまうことです。

使いどころとトレードオフ

MCMCの強みは汎用性です。解析解のない複雑なモデル — 階層モデル、混合分布、非線形モデル — でも、事後密度に比例する量が書ければ推論できます。代価は計算時間で、1回の推論に数分から数時間かかることも珍しくありません。速度が最優先の場面では、事後分布を単純な分布で近似する変分推論という代替もありますが、近似の粗さ(とくに不確実性を小さく見積もりがちなこと)と引き換えです。「正確だが遅いMCMC、速いが近似の変分推論」という使い分けが現在の標準的な構図で、モデル開発時はMCMCで確かめ、大規模運用では近似手法に切り替える、といった組み合わせも実務ではよく行われます。

試問

等高線で描いた「二つの山」の分布から、メトロポリス法の歩き手で標本を集めます。歩き手は毎歩、いまの位置の近くにランダムな候補を出し、山の高いほうへはほぼ必ず移り、低いほうへはときどきだけ移ります(候補へ移れた割合が受理率)。はじめの200歩は出発点の影響が残るので捨て(burn-in)、その後の足跡を標本として数えます。候補を出す範囲の広さ=歩幅は数値だけの選択肢で、いちばん小さい歩幅から始まります。

✦ 命題挑戦中

歩幅を選んで連鎖を走らせ、二つの山を目標どおりの割合で行き来する標本を集めよ

  • 目標「二つの山」で、burn-in 後の標本を 3000 歩ぶん以上集める(いま 0 歩)
  • 受理率を 15〜60% に保つ(いま —)
  • 右の山(x>0)にいた割合が目標 45% の ±15 ポイント以内(いま —)
⚡ 体験: メトロポリス法で分布の上を歩く
目標
歩幅
目標分布と歩き手0 歩トレース(x 座標)burn-in-4-2024歩数 →標本
歩数0受理率—x>0 にいた割合—目標45%

金の歩き手は分布のすそ野の外れから出発します(つかんで置き直せます)。「1歩」を押すと、近くに候補を出し、密度の比 r で受理・棄却を決めます。

等高線 = 目標分布(正規化定数は知らないまま)。金 = 受理された足跡、朱 = 棄却された候補、白い点の雲 = 集めた標本