ベイズ校正 MCMC シミュレーター

Q: 有効サンプル数 ESS が小さいとなぜ問題ですか？

MCMC のサンプルは互いに相関しているため、5000 サンプル取っても「独立サンプル換算」では数百〜数十ということがよくあります。これが ESS です。事後平均の標準誤差は σ_post/√ESS で決まるため、ESS が 100 を切ると推定値の標準誤差が事後 std の 10% 以上となり、結果の小数第2位以下が信頼できません。ESS を増やすには (1) 提案幅を最適化する、(2) サンプル数を増やす、(3) Hamiltonian Monte Carlo や NUTS など効率的なサンプラに切り替える、の三つが選択肢です。

ベイズ校正 MCMC シミュレーター

事前分布・尤度・観測データを与えて、CAEモデルのパラメータをベイズ流に校正するツールです。共役正規モデルの解析事後分布をリアルタイム計算し、Metropolis–Hastings法の受理率と有効サンプル数 (ESS) を同時に評価して、UQ／モデル較正の感覚をつかめます。

パラメータ設定

事前平均 μ_p

パラメータ θ の事前知識の中心

事前標準偏差 σ_p

大きいほど弱情報的（データに従う）

観測雑音 σ_y

観測モデルの既知標準偏差

観測平均 ȳ

N 個の観測の標本平均

観測数 N

事後 std は 1/√N で縮む

MH 提案幅 step

2.38·σ_post 付近が最適

MCMC サンプル数

バーンイン

初期遷移サンプルの破棄数

計算結果

—

事後平均 μ_post

—

事後標準偏差 σ_post

—

95% 信用区間幅

—

MH 受理率 (%)

—

有効サンプル数 ESS

—

情報利得 (nats)

—

事前・尤度・事後重ね描き + MH walker

青：事前 p(θ)、赤：尤度 p(y|θ)、緑：事後 p(θ|y)。白い丸は事後分布上をランダムウォークする MH チェーンの現在位置です。

事前 / 尤度 / 事後 PDF

MCMC trace plot — θ vs iteration

理論・主要公式

$$p(\theta|y) \propto p(y|\theta)\,p(\theta),\qquad \alpha = \min\left(1,\, \frac{\pi(\theta^{*})}{\pi(\theta_t)}\right)$$

π は無正規化の事後密度。提案 θ* を確率 α で受理し、十分長い chain で π からのサンプルが得られる（Metropolis–Hastings）。

$$\frac{1}{\sigma_{\text{post}}^{2}} = \frac{1}{\sigma_p^{2}} + \frac{N}{\sigma_y^{2}},\qquad \mu_{\text{post}} = \sigma_{\text{post}}^{2}\!\left(\frac{\mu_p}{\sigma_p^{2}} + \frac{N\,\bar{y}}{\sigma_y^{2}}\right)$$

正規－正規共役の解析事後分布。観測数 N が増えると 1/σ_post² が線形に増え、σ_post は 1/√N で縮む。

$$\text{CI}_{95} = \mu_{\text{post}} \pm 1.96\,\sigma_{\text{post}},\qquad I = \log\!\frac{\sigma_p}{\sigma_{\text{post}}}\ [\text{nats}]$$

95% 信用区間と情報利得 I（事前→事後で減少した不確かさを nats で表した量）。

ベイズ校正（MCMC）— Metropolis–Hastings

🙋

ベイズ校正って、CAE のモデルを実験データに合わせる作業ですよね？普通の最小二乗で合わせるのとどう違うんですか？

🎓

いい質問だね。最小二乗（最尤）は「データに一番合う1点の θ」を返すだけだけど、ベイズはパラメータ θ を確率分布として扱うのがポイント。事前分布 p(θ) に過去の知見や妥当な範囲を入れて、観測の尤度 p(y|θ) と掛け算して事後分布 p(θ|y) を作る。結果は「点」じゃなく「ばらつき付きの分布」になるから、設計余裕や信頼区間を自然に語れるんだ。例えば材料のヤング率を 200 ± 5 GPa と推定できれば、後段の応力計算に幅で持ち込める。

🙋

なるほど。それで、なぜわざわざ MCMC が出てくるんですか？事後分布って解析的に書けないんでしょうか？

🎓

このツールみたいに事前も尤度も正規だと「共役」と呼ばれる特殊なケースで、事後分布が閉形式で出る。だから解析事後（緑のカーブ）と MCMC 推定が一致して気持ちいい。でも実際の CAE モデル — 例えば乱流の k-ε パラメータや非線形材料の硬化則 — では尤度 p(y|θ) を評価するのにフルの数値解析が要る。閉形式は無理だから、π(θ) を「比 π(θ*)/π(θ_t)」の形で計算しながら受理／棄却で進む MH が万能ツールとして使われるんだ。

🙋

MH 提案幅 step を動かすと、受理率がギュッと変わりますね。最適な値ってどう決めるんですか？

🎓

古典的な結果として、Roberts–Gelman–Gilks (1997) が「1次元なら約 44%、高次元なら約 23.4%」が最適と示している。経験則は σ_proposal ≒ 2.38·σ_post。step を大きくしすぎると受理率がほぼ 0% に落ちて鎖がフリーズし、小さくしすぎると受理率は高くても自己相関が強くて有効サンプル数 ESS が減る。受理率 30〜50% を狙ってチューニングするのが実務感覚だね。デフォルトの 0.15 は σ_post ≒ 0.07 に対して比 ≈ 2.1 なので、ほぼ最適点に当ててある。

🙋

最後に「情報利得」って何ですか？単位が nats って初めて見ました。

🎓

事前分布と事後分布の差を、情報量の単位で測ったものだよ。今回は近似として log(σ_p/σ_post) を出してるけど、これは「観測でどれだけ不確かさが縮んだか」の指標。底を e にすると単位は nats（底 2 なら bits）。デフォルトでは log(2/0.0707) ≒ 3.34 nats、bits 換算で約 4.8 bits。つまり「観測で約 5 ビット分、θ について新しい情報が得られた」ということ。CAE の校正実験の「情報的価値」を比較するときに便利だよ。

よくある質問

最小二乗法（最尤推定）はパラメータを1点の値として推定しますが、ベイズ校正はパラメータを確率分布として扱い、事前知識 p(θ) と観測 p(y|θ) を組み合わせて事後分布 p(θ|y) を出します。これにより「最適値」だけでなく「不確かさの幅」と「他の値である確率」まで定量化できます。CAEモデルの校正では観測データが少ない・ノイズが多いケースが多く、点推定だけでは過信のもとです。本ツールでは観測数 N を 5〜500 の範囲で動かすと、事後標準偏差 σ_post が 1/√N で縮む様子が直接見えます。

Roberts–Gelman–Gilks (1997) の漸近最適化に基づき、1次元では受理率 ≈ 44%、高次元では ≈ 23.4% が最適とされます。経験則としては提案標準偏差を事後標準偏差の約 2.38 倍にすると受理率が最適値付近に落ち着きます。本ツールではこの目安からの乖離を bell-shape で受理率に反映しており、step を大きくしすぎると受理率がほぼ 0% に落ちて鎖が動かなくなり、小さくしすぎると受理率は高くなるものの相関が強く有効サンプル数 ESS が減ります。step を動かしながら受理率 30〜50% の範囲を探すのが実務的なチューニングです。

MCMC の初期値は事後分布から外れていることが多く、最初の数百〜数千ステップは事後の高密度域に「移動中」の遷移サンプルが続きます。これを推定に混ぜると事後平均や分散にバイアスが乗るため、初期のサンプルを破棄します。これがバーンインです。本ツールのデフォルト 500 ステップは弱情報的な事前分布と数千サンプルの鎖に対して保守的な値です。trace plot を見て「左端のドリフトが消える点」を目視で確認すると安全です。さらに収束を厳密に判定したい場合は複数チェーンを走らせ Gelman–Rubin の R̂ < 1.01 を確認します。

MCMC のサンプルは互いに相関しているため、5000 サンプル取っても「独立サンプル換算」では数百〜数十ということがよくあります。これが ESS です。事後平均の標準誤差は σ_post/√ESS で決まるため、ESS が 100 を切ると推定値の標準誤差が事後 std の 10% 以上となり、結果の小数第2位以下が信頼できません。ESS を増やすには (1) 提案幅を最適化する、(2) サンプル数を増やす、(3) Hamiltonian Monte Carlo や NUTS など効率的なサンプラに切り替える、の三つが選択肢です。

実世界での応用

CFD／FEM コードのパラメータ校正：乱流モデルの定数（k-ε の C_μ、Cε1、Cε2 など）、非線形材料の硬化則パラメータ、コンクリート損傷モデルの限界ひずみなど、解析者が「えいや」で決めていた定数を、実験データから事後分布で推定する用途です。NASA や DOE は不確かさ定量化（VVUQ：Verification, Validation, Uncertainty Quantification）の必須プロセスとして MCMC ベイズ校正を採用しています。本ツールの正規−正規モデルは、その入門として最もシンプルな共役ケースです。

気候・地球科学モデルの校正：気候感度や雲フィードバック係数のように直接観測できないパラメータを、過去の気温トレンドや衛星観測から推定する場面で MCMC が標準的に用いられます。Hadley Centre や NOAA の解析報告書では事後分布のヒストグラムが必ず示されます。観測が少なく事前分布の影響が大きい領域では、本ツールで σ_p を動かしたときに事後がどれだけ事前に引っ張られるかが直感的に分かります。

機械学習ベイズ推論：ベイズニューラルネット（BNN）、ガウス過程の超パラメータ推定、A/B テストの効果量推定など、ML の現場でも MCMC は活躍します。最近は HMC／NUTS（Stan、PyMC、NumPyro）が主流ですが、その基礎には本ツールが扱う Metropolis–Hastings の受理／棄却の発想があります。受理率と ESS を見る「鎖の健康診断」は HMC でも同じです。

医薬品・薬物動態（PK/PD）モデル：体内動態のコンパートメントモデルに対し、少数の被験者データから個人差込みの事後分布を推定する場面で、ベイズ階層モデル＋MCMC が業界標準です。FDA の審査資料にも事後分布の信用区間が必須項目として並びます。本ツールの「観測数 N が小さいと σ_post が大きい」現象は、被験者数の少ない第I相試験そのものの構造を表しています。

よくある誤解と注意点

まず最大の落とし穴が、「受理率が高い＝良いサンプラ」だと思い込むことです。提案幅 step を極端に小さくすれば受理率は 90% を超えますが、隣のサンプルとほぼ同じ場所にしか動かないため、5000 サンプル取っても独立サンプル換算 ESS は数十、ということが起こります。本ツールで step=0.05 を試してから ESS の小ささを確認してみてください。受理率は単独で見ず、必ず ESS とセットで評価することが重要です。

次に、「事前分布を平均0・分散大の弱情報にしておけば常に客観的」という誤解。事前を広くしすぎると、観測が少ないときに事後が「データに引っ張られすぎ」て不安定になります（典型的にはばらつきや極端値が出やすい）。逆に、強情報的な事前（σ_p を小さく）はモデルへの自信を表すと同時に、観測との不一致を隠してしまうリスクがあります。本ツールで N=5、σ_p=0.5 と N=5、σ_p=10 を比較すると、事前の影響度が体感できます。事前は「決め打ち」ではなく、複数のシナリオを試して感度を見るのがプロの作法です。

最後に、「MCMC が収束したかどうかを 1 本の鎖の見た目だけで判断する」こと。1 本の鎖が一見「安定」に見えても、別の初期値から走らせると全く別のモードに収束していることがあります（特に多峰性事後分布）。実務では複数チェーンを並走させ、Gelman–Rubin の R̂ がすべて 1.01 未満であることを確認します。本ツールは 1 次元 1 鎖の概念学習用で、R̂ プロキシを参考表示します。実プロジェクトでは PyMC や Stan の summary() で R̂・ESS・MCSE を必ず確認してください。

ベイズ校正（MCMC）— Metropolis–Hastings

よくある質問

実世界での応用

よくある誤解と注意点

使い方ガイド

具体的な計算例

実務での注意点