パラメータ設定
乱数は線形合同法（LCG）で seed=42 固定。同じ設定なら毎回同じ結果になります。
サンプルチェーン・ヒストグラム・自己相関
上＝トレース x_t（青）／中＝ヒストグラム（緑）と真の密度（黒）／下＝自己相関 ρ_k（ラグ 1〜50）
理論・主要公式
メトロポリス・ヘイスティングス法は、ガウス提案分布 $q$ から候補 $x'$ を生成し、受理確率 $\alpha$ で受理／棄却します。提案 $q$ が対称なら $q(x|x') = q(x'|x)$ がキャンセルされ、確率密度の比だけで判定できます。
提案：
$$x' \sim \mathcal{N}(x_t, \sigma_\text{prop}^2)$$
受理確率：
$$\alpha = \min\!\left(1, \frac{\tilde p(x')}{\tilde p(x_t)}\right)$$
更新規則（$u \sim U(0,1)$）：
$$x_{t+1} = \begin{cases} x' & (u < \alpha) \\ x_t & (\text{otherwise}) \end{cases}$$
$\tilde p(x)$ は正規化されていない密度でかまいません（比だけ使うため）。1次元ガウス目標では受理率が約 44％ になる $\sigma_\text{prop}$ が最適です。
MCMCメトロポリス法サンプラーとは
🙋
確率分布からサンプリングしたいだけなのに、なぜマルコフ連鎖なんていう難しそうな方法を使うんですか？普通に乱数から作れないんですか？
🎓
ざっくり言うと、複雑な分布、特に正規化定数（積分が1になるようにする定数）が計算できない分布からは、直接サンプリングできないからなんだ。たとえばベイズ推定の事後分布は分母の周辺尤度が高次元の積分で計算困難なことが多い。そこでメトロポリス・ヘイスティングス法は「比だけ」で動く——$p(x')/p(x)$ の比は正規化定数がキャンセルされるから、密度の形だけ分かれば良い。シミュレーターの「目標分布タイプ」を1（2峰混合）にして、ヒストグラムが2つの山を再現していくのを見てみて。
🙋
提案分散 σ_prop を 0.1 にしたら、トレースがまっすぐな線みたいになって、ヒストグラムが片方の山しか描けてないんですけど。これってバグですか？
🎓
バグじゃなくて、これがMCMCの「混合（mixing）」問題の典型例だよ。σ_prop が小さすぎると一歩ずつしか動けないから、谷を越えてもう一つの山にジャンプするのに何千ステップもかかる。受理率カードを見ると 95％ 以上になってるはずだけど、これは「ほぼ毎回受理してるけど、ほとんど動いてない」状態。逆に σ_prop=5 にすると、ほとんど棄却されてチェーンが同じ場所に張り付く。受理率が 30〜50％ くらいになる σ_prop が「ちょうどいい」目安なんだ。
🙋
なるほど！じゃあ最適な σ_prop さえ選べば完璧にサンプリングできるってことですか？
🎓
実はそこにもう一つ罠があって、隣り合うサンプルは似てるから、N=2000 個取っても情報量は独立 N個分にならない。これを表すのが「有効サンプル数（ESS）」で、自己相関の下のグラフのバーが速く0に近づくほど ESS が大きくなる。CAEの不確実性定量化では、目標精度に必要なESSを先に決めて、それから連鎖長を決めるのが普通だよ。試しに σ_prop を最適値から外して、ESS のカードがどれくらい減るか見てみると感覚がつかめると思う。
よくある質問
マルコフ連鎖の定常分布が目標 p(x) になるための「詳細釣り合い条件 p(x)q(x'|x)α(x,x') = p(x')q(x|x')α(x',x)」を満たす最大の受理率がこの形です。提案 q が対称なら q(x'|x)=q(x|x') なので、比 p(x')/p(x) だけで決まります。これがメトロポリス法の核心で、正規化定数が不明な分布でもサンプリングできる魔法の源です。
Roberts, Gelman, Gilks (1997) によるガウス目標の漸近解析で、有効サンプル数を最大化する最適受理率が1次元では 0.44、多次元では次元数が大きくなるにつれ 0.234 に収束することが示されました。これより受理率が高すぎると一歩が小さすぎて混合が遅く、低すぎると棄却が多くチェーンが止まります。実務では 0.2〜0.5 を目安に σ_prop を調整します。
山と山の間の谷の確率密度が低いと、提案点がそこに落ちても確率比が極端に小さく棄却されるため、もう一方の山へ「ジャンプ」できなくなります。σ_prop を大きくして谷を越えやすくする、複数の初期値で並列に走らせる、温度を上げて谷を浅くするレプリカ交換法（パラレル・テンパリング）を使うなどの対策があります。
トレースプロットを見て、初期値の影響が消え水平に揺らぐ領域に入った点までを目安にします。本ツールではデフォルトで 200 ステップに設定していますが、複雑な分布や悪い初期値ではこれより長く必要なこともあります。診断指標としては Gelman-Rubin の R̂ (複数チェーンの分散比) や、自己相関時間の数倍程度を取るのが定番です。
実世界での応用
ベイズ推定とパラメータ推定：事後分布 $p(\theta|D) \propto p(D|\theta)p(\theta)$ から $\theta$ をサンプリングする最も標準的な手法です。観測データから物理モデルのパラメータ（例：材料の弾性係数の事後分布）を推定する場面で、解析的に積分できない事後分布も MCMC なら扱えます。Stan や PyMC など主要なベイズ推定ライブラリの中核アルゴリズムです。
統計力学とイジングモデル：古典的にメトロポリス法はモンテカルロ統計力学で、エネルギー $E$ を持つ状態のボルツマン分布 $p \propto e^{-E/kT}$ からサンプリングするために発明されました（Metropolis et al., 1953）。スピン系の相転移、合金の規則-不規則転移、タンパク質フォールディングなど、平衡状態の物理量を計算するのに使われ続けています。
CAEの不確実性定量化（UQ）：有限要素解析の入力パラメータ（材料定数、境界条件、幾何公差）が確率分布を持つとき、出力（応力、変形）の事後分布や信頼性指標を MCMC で推定します。FORM/SORM では捉えきれない複雑な失敗領域や、複数のモードを持つ分布の評価に有効です。サロゲートモデル（クリギング、Gaussian Process）と組み合わせて計算コストを下げるのが現代的なやり方です。
機械学習と確率的プログラミング：潜在変数モデル、混合モデル、階層モデルなど、最尤推定では局所解にはまる複雑なモデルの学習に MCMC が使われます。トピックモデル（LDA）の Collapsed Gibbs サンプラーや、ベイジアンニューラルネットの重みのサンプリング、強化学習の方策評価などにも応用されています。
よくある誤解と注意点
最も多い誤解は、「受理率が高いほど良い」と思ってしまうことです。シミュレーターで σ_prop=0.1 にすると受理率は 95％ を超えますが、これは「ほぼ動いていない」だけで、サンプルの質はむしろ最悪です。逆に σ_prop=5 にすると受理率は 10％ 程度に下がりますが、これも頻繁に棄却されてチェーンが固まる悪い状態です。1次元では受理率 40〜50％、多次元では 20〜30％ を目指して σ_prop を調整するのが正解で、本ツールではデフォルトの σ_prop=1.0 がほぼ最適値になっています。
次に、サンプル数 N と独立サンプル数を混同する誤りです。MCMCのサンプル $x_1, x_2, \ldots, x_N$ は前のサンプルに依存するため、独立ではありません。たとえば自己相関時間が $\tau=20$ なら、N=2000 のチェーンの実効的な情報量はわずか 100 個程度です。本ツールの ESS カードはこれを簡易計算で示しています。論文で「N=10000のMCMCを実行」と書くだけでは不十分で、必ず ESS や $\hat R$ などの診断指標を併記すべきです。
最後に、バーンインを切り捨てれば収束保証が得られると思い込む誤解です。バーンインは「初期値の影響を消す」ためのもので、「定常分布に達したことの保証」ではありません。複数の山がある分布で 1 つの山にしか到達していない場合、いくらバーンインを長くしても他の山は見えません。実務では複数の初期値から並列にチェーンを走らせ、Gelman-Rubin 統計量 $\hat R$ が 1 に近いことを確認する、トレースプロットを目視で点検する、といった複数の診断を組み合わせて収束を判断します。