ベイズキャリブレーション
理論と物理
ベイズキャリブレーションとは
先生、ベイズキャリブレーションって、実験データでモデルのパラメータを推定するんですか？ 普通の最適化と何が違うんですか？
いい質問だね。ざっくり言うと、パラメータを「1つの値」ではなく「確率分布」として推定するのがベイズキャリブレーションだ。通常の最小二乗法や最適化は「実験データに一番合う値はこれだ！」と点推定するけど、ベイズは「この範囲にこのくらいの確率でいる」という分布で答えを返す。
分布で答えが返ってくると、何が嬉しいんですか？
たとえば自動車の衝突解析でFEMモデルのヤング率を推定するとしよう。最適化だと $E = 210\,\text{GPa}$ という1つの値しか出ない。でもベイズなら「$E$ は平均 $210\,\text{GPa}$ で標準偏差 $5\,\text{GPa}$ の分布」と出る。この不確かさ情報を後段の信頼性解析にそのまま渡せるんだ。
さらに大事なのは、エンジニアの経験知（事前分布）を正式に組み込めること。「この材料のヤング率はだいたい $200 \sim 220\,\text{GPa}$ だろう」という事前知識と、実験データの両方を活かしてパラメータを絞り込む。これがベイズの最大の強みだよ。
経験知を数式に取り込めるって、すごく実務的ですね！ でもそれって具体的にはどういう数学なんですか？
ベイズの定理と基本方程式
ベイズキャリブレーションの数学的骨格はベイズの定理だ。パラメータベクトル $\boldsymbol{\theta}$ と観測データ $\mathbf{D}$ に対して：
各項の物理的意味
- $p(\boldsymbol{\theta}|\mathbf{D})$（事後分布）：実験データを見た後のパラメータの確率分布。これがキャリブレーションの最終出力。
- $p(\mathbf{D}|\boldsymbol{\theta})$（尤度関数）：あるパラメータ値のもとで実験データが得られる確率。モデル予測と実験データの「ズレ具合」を測る。
- $p(\boldsymbol{\theta})$（事前分布）：データを見る前のパラメータに関する知識。文献値、経験、物理的制約を反映する。
- $p(\mathbf{D})$（エビデンス/周辺尤度）：正規化定数。モデル比較に使うが、パラメータ推定では定数として無視できる。
尤度関数って、具体的にはどう書くんですか？
一番シンプルなケースを考えよう。$N$ 個の実験データ $d_i$ に対して、モデル予測 $M(\mathbf{x}_i, \boldsymbol{\theta})$ との差が独立同分布の正規誤差に従うと仮定すると：
ここで $\sigma$ は観測ノイズの標準偏差だ。これを既知として固定してもいいし、$\boldsymbol{\theta}$ と一緒にベイズで推定してもいい。実務では $\sigma$ も推定対象に含める方が堅実だよ。誤差の大きさが正確にわからないケースが多いからね。
でも先生、FEMモデルって物理をある程度近似してますよね？ モデルと実験が「正規誤差だけ」で説明できるって仮定は甘くないですか？
鋭いね！ まさにそこを体系的に扱うのがKennedy-O'Haganフレームワークだ。
Kennedy-O'Haganフレームワーク
2001年にKennedyとO'Haganが提案した枠組みは、CAEのベイズキャリブレーションにおける事実上の標準になっている。核心は、実験値を次のように分解することだ：
Kennedy-O'Hagan各項の意味
- $M(\mathbf{x}, \boldsymbol{\theta})$（シミュレーション出力）：入力条件 $\mathbf{x}$ とキャリブレーションパラメータ $\boldsymbol{\theta}$ を与えたときのFEM/CFD等の計算結果。
- $\delta(\mathbf{x})$（モデル不整合項）：パラメータをどう調整しても埋められないモデル自体の構造的な誤差。メッシュ離散化誤差、省略された物理現象、理想化された境界条件などに起因する。通常、ガウス過程（GP）で表現する。
- $\varepsilon$（観測誤差）：実験計測に伴うランダムノイズ。$\varepsilon \sim \mathcal{N}(0, \sigma_\varepsilon^2)$ と仮定するのが一般的。
$\delta(\mathbf{x})$ を入れるのと入れないのでは、何がどう変わるんですか？
$\delta$ を入れないと、モデルの構造的誤差を全部パラメータ $\boldsymbol{\theta}$ に押し付けることになる。たとえば、実際には境界条件の不備が原因でモデルと実験がずれているのに、ヤング率を極端な値にして辻褄を合わせてしまう。これがパラメータ補償と呼ばれる問題で、物理的に意味のないパラメータ値が推定される原因だ。
$\delta(\mathbf{x})$ を明示的に入れることで、「モデルでは再現できない部分」を分離して、パラメータの推定を物理的に妥当な範囲に保てるんだ。
なるほど…！ ちゃんとモデルの限界を認めた上で推定するってことですね。でも $\delta$ 自体はどうやって決めるんですか？
$\delta(\mathbf{x})$ はガウス過程（GP）としてデータから学習する。つまり、$\delta$ の形状は事前に決めるのではなく、実験データとモデル予測の残差パターンからベイズ推論で自動的に推定されるんだ。GPのカーネル関数（たとえば二乗指数カーネル）のハイパーパラメータも一緒に推定する。
事前分布の選び方
事前分布ってどうやって決めるんですか？ 主観的で恣意的になりませんか？
事前分布の設定は確かにベイズの最も議論が分かれるポイントだ。実務では次のような使い分けが多い：
- 情報的事前分布：材料データベースの値、過去の実験結果、ハンドブック値をもとに設定。例：構造用鋼のヤング率なら $E \sim \mathcal{N}(210, 5^2)\,[\text{GPa}]$。
- 弱情報的事前分布：物理的に妥当な範囲を緩く制約。例：摩擦係数は $\mu \sim \text{Uniform}(0.05, 0.8)$。
- 無情報的事前分布：Jeffreysの事前分布など、データに語らせたいとき。ただし高次元では事後分布の収束が遅くなる。
現場で多いのは弱情報的事前分布だね。「物理的にありえない領域は排除するが、データに十分な影響力を持たせる」というバランスだ。
事前分布の選び方で結果が大きく変わることもあるんですか？
データが少ないときは事前分布の影響が大きく出る。逆にデータが十分にあれば、事前分布が多少違っても事後分布はほぼ同じところに収束する。これをベイズの漸近的一致性という。だから実務上は、事前分布の感度分析（事前分布を変えて結果がどの程度変わるかを確認）を必ずやるべきだよ。
数値解法と実装
MCMC法の基礎
ベイズの定理自体はシンプルですけど、実際に事後分布を計算するのは大変そうですよね？ パラメータが多いと積分できないと思うんですが…
その通り。事後分布の正規化定数 $p(\mathbf{D})$ は高次元積分なので解析的に求められない。そこで使うのがマルコフ連鎖モンテカルロ（MCMC）法だ。事後分布から直接サンプルを生成する手法で、サンプル数を増やせば事後分布を任意の精度で近似できる。
MCMCってどういう仕組みなんですか？
最も基本的なMetropolis-Hastingsアルゴリズムの流れはこうだ：
- 初期値 $\boldsymbol{\theta}^{(0)}$ を設定
- 提案分布 $q(\boldsymbol{\theta}^*|\boldsymbol{\theta}^{(t)})$ から候補点 $\boldsymbol{\theta}^*$ を生成
- 採択確率を計算：
確率 $\alpha$ で候補 $\boldsymbol{\theta}^*$ を採択し、$(1-\alpha)$ で棄却して現在の値に留まる。これを何千〜何万回も繰り返すと、サンプル列が事後分布に収束するんだ。
イメージとしては、山の中をランダムウォークしながら、標高の高い（尤度×事前確率が大きい）場所に長く滞在する仕組みだよ。
おお、なるほど！ でもFEMの計算って1回に数分〜数時間かかりますよね？ 何万回も回すのは現実的なんですか？
鋭い。まさにそこがCAEでベイズキャリブレーションを使う最大のボトルネックだ。対策は後で詳しく話すサロゲートモデルだけど、まずは各MCMCアルゴリズムの特徴を整理しよう。
MCMCアルゴリズムの比較
| アルゴリズム | 特徴 | 適用場面 | 典型的な採択率 |
|---|---|---|---|
| Metropolis-Hastings | 最も基本的。実装が簡単 | 低〜中次元（〜10パラメータ） | 20〜50% |
| Adaptive Metropolis (AM) | 共分散行列を自動調整 | 中次元。チューニング不要 | 23%（最適） |
| DRAM | Delayed Rejection + AM。棄却時に再提案 | 中次元。多峰性に一定の対応力 | 30〜50% |
| ハミルトニアンMC (HMC) | 勾配情報を利用。高次元でも効率的 | 高次元（100+）。微分可能モデル | 65〜90% |
| NUTS | HMCの自動チューニング版。Stan標準 | 高次元。ハンズフリー | 80〜95% |
| Ensemble Sampler (emcee) | 多数のウォーカーが並列探索 | 中次元。並列化容易 | 20〜50% |
| TMCMC / SMC | 逐次モンテカルロ。多峰性に強い | 多峰性事後分布。モデル比較 | — |
こんなにたくさんあるんですね！ CAEでよく使われるのはどれですか？
Sandia国立研究所のDakotaではDRAMが標準実装されていて、実務で最もよく使われている。勾配情報が取れるならHMC/NUTSが圧倒的に効率いい。Pythonで手軽にやるならPyMC（NUTSベース）かemceeだね。
収束診断
MCMCが「収束した」っていうのはどうやって判断するんですか？
これは非常に重要なポイントだ。収束してないMCMCの結果を使うと、全く意味のない事後分布が出てくる。標準的な収束診断は以下の通り：
- $\hat{R}$（Gelman-Rubin統計量）：複数の独立チェーンの分散比。$\hat{R} < 1.01$（厳密には $< 1.05$）なら収束と判断。
- トレースプロット：サンプル列の時系列をプロット。「毛虫のような」混合の良いパターンが理想。トレンドやジャンプがあれば未収束。
- 有効サンプルサイズ (ESS)：自己相関を考慮した実質的なサンプル数。ESS $> 400$（各パラメータで）が目安。
- バーンイン除去：初期の過渡的なサンプル（通常チェーン長の10〜50%）を捨てる。
$\hat{R} < 1.01$ を確認して、トレースプロットが安定していて、ESSが十分ならOKということですね。。
サロゲートモデルによる高速化
さっき「FEMの計算を何万回回すのは大変」って話でしたけど、実際どうやって解決するんですか？
最も一般的な解決策がサロゲートモデル（代理モデル）だ。FEMの入出力関係を高速な近似モデルで置き換える。主な手法は：
- ガウス過程回帰（GP/Kriging）：Kennedy-O'Haganフレームワークとの親和性が高い。予測の不確かさも出る。
- 多項式カオス展開（PCE）：入力パラメータの確率分布を直交多項式で展開。不確かさ伝播と同時に使える。
- ニューラルネットワーク：高次元・非線形が強いケースで有効。ただしデータ量が必要。
たとえばFEMを100〜500回実行してGPを学習させれば、MCMCの各ステップではGP予測（ミリ秒オーダー）だけで済む。原理的にはFEMの計算回数を $10^4 \to 10^2$ に削減できるんだ。
サロゲートの精度が悪いと、キャリブレーション結果も間違いになりませんか？
その通り。だからサロゲートの検証が必須だ。具体的には、Leave-one-out交差検証で予測精度を確認する。 $Q^2 > 0.99$ が目安。もう一つの戦略として、逐次サロゲート更新がある。MCMCの途中で事後分布の高確率領域にFEM計算点を追加して、サロゲートを精密化していく。これで効率と精度の両立ができる。
実践ガイド
キャリブレーションのワークフロー
実際にベイズキャリブレーションをやるときの手順を教えてください！
典型的なワークフローは以下の7ステップだ：
- キャリブレーション対象パラメータの選定：感度分析（Morris法やSobol指数）で影響度の大きいパラメータを絞る。パラメータ数は5〜10以下に抑えたい。
- 実験データの準備：測定値、測定条件、測定誤差の情報を整理。データの品質が結果を左右する。
- 事前分布の設定：各パラメータの範囲と分布形状を決定。物理的制約（非負、上限など）を反映。
- 尤度関数の構築：誤差モデル（正規、対数正規、t分布など）を選択。Kennedy-O'Haganの $\delta$ を入れるか判断。
- サロゲートモデルの構築：実験計画法（ラテン超方格など）でFEMを実行し、GPやPCEを学習。
- MCMCの実行と収束診断：複数チェーンを並列実行。$\hat{R}$、ESS、トレースプロットで確認。
- 事後分布の解析：周辺分布、相関、予測区間を可視化。事後予測が実験データと整合するか検証。
実践例：溶接残留応力モデル
具体的な例を見たいです！ どんなケースでベイズキャリブレーションが活きるんですか？
典型的な例として、溶接の熱弾塑性FEM解析を考えよう。残留応力を精度よく予測したいが、いくつかのパラメータに大きな不確かさがある：
- 溶接入熱効率 $\eta$：文献値は $0.6 \sim 0.9$ と幅がある
- 熱源形状パラメータ（Goldakモデルの $a_f, a_r, b, c$）
- 高温域の降伏応力（実測困難）
これらをキャリブレーション対象とし、中性子回折やX線回折で測った残留応力データ（たとえば溶接線直交方向の10点）を実験値として使う。
溶接の入熱効率を「0.75」ってピンポイントで決める代わりに、「0.6〜0.9の範囲で実験データと照らし合わせて確率分布を求める」ってことですね？
そうだ。そしてキャリブレーション後に「$\eta$ は平均 $0.78$、95%信用区間 $[0.72, 0.84]$」みたいな結果が得られる。さらに事後分布を使って残留応力の予測区間も出せるから、品質保証や規格適合の判断に使えるんだ。
よくある失敗と対策
ベイズキャリブレーションで「ハマりがち」なポイントってありますか？
現場で本当によく見る失敗パターンをまとめるよ：
| 失敗パターン | 症状 | 対策 |
|---|---|---|
| パラメータの非識別性 | 事後分布が事前分布とほとんど変わらない。または複数パラメータが強く相関 | 感度分析で影響度の低いパラメータを固定する。実験の種類を増やす |
| モデル不整合の無視 | パラメータが物理的にありえない値に推定される | Kennedy-O'Haganの $\delta(\mathbf{x})$ を導入 |
| MCMCの未収束 | $\hat{R} > 1.1$、トレースプロットが漂流 | チェーン数・ステップ数を増やす。提案分布のスケールを調整 |
| サロゲートの精度不足 | 事後分布にアーティファクト（偽のピーク等） | 事後高確率領域に計算点を追加してサロゲートを再学習 |
| 尤度関数の誤設定 | 事後予測の信用区間が実験データをカバーしない | 誤差モデルを見直す。$\sigma$ を推定対象にする |
| 事前分布が狭すぎる | 真値が事前分布の裾にあり、事後分布が裾に張り付く | 事前分布を広げる。事前分布感度分析を行う |
パラメータの非識別性ってよくあるんですか？
めちゃくちゃ多い。たとえば熱伝導問題で熱伝導率 $k$ と対流熱伝達係数 $h$ を同時にキャリブレーションしようとすると、温度データだけでは2つを分離できないことがある。結果として $k$ と $h$ の間に強い負の相関が出て、「$k$ が大きくて $h$ が小さい」でも「$k$ が小さくて $h$ が大きい」でもデータを説明できてしまう。こういうときは、一方を文献値で固定するか、異なる種類の実験データ（表面温度＋内部温度など）を追加するしかない。
ソフトウェア比較
ツール一覧と機能比較
ベイズキャリブレーションを実行するには、どんなソフトウェアが使えるんですか？
主要ツールをまとめよう。大きく分けて「UQ専用フレームワーク」と「汎用ベイズ推論ライブラリ」の2種類がある：
| ツール | 開発元 | MCMC手法 | サロゲート | Kennedy-O'Hagan | ライセンス |
|---|---|---|---|---|---|
| Dakota | Sandia国立研究所 | DRAM, MH, DREAM | GP, PCE, NN | 対応 | OSS (LGPL) |
| UQLab | ETH Zurich | AIES, AM, HMC | GP, PCE, SVR | 対応 | 教育無料 / 商用有料 |
| MUQ | MIT | MH, DILI, SMC | GP | 対応 | OSS (BSD) |
| PyMC | コミュニティ | NUTS, MH, SMC | 外部連携 | 手動実装 | OSS (Apache) |
| Stan | コミュニティ | NUTS (HMC) | 外部連携 | 手動実装 | OSS (BSD) |
| emcee | D. Foreman-Mackey | Affine Invariant | 外部連携 | 手動実装 | OSS (MIT) |
| Ansys optiSLang | Ansys Inc. | 独自実装 | MOP, GP | 限定的 | 商用 |
| SIMULIA Isight | Dassault Systemes | 独自実装 | RSM, Kriging | 限定的 | 商用 |
選定の指針
結局、どれを使えばいいんですか？
使い分けの目安はこんな感じだ：
- CAEソルバーとの統合重視 → Dakota。AbaqusやOpenFOAMとのインターフェースが充実。
- MATLABユーザー → UQLab。UQ全般が一通りできて、ドキュメントも充実。
- Pythonで柔軟にやりたい → PyMC + GPy（サロゲート）。論文のプロトタイプには最適。
- 高次元パラメータ → Stan（NUTS）。勾配ベースで100+パラメータもいける。
- 既存の商用ワークフローに組み込みたい → optiSLangやIsight。GUIベースで非専門家も使いやすい。
まずはDakotaかPyMCで始めてみるのがよさそうですね！
先端技術
最新の研究動向
ベイズキャリブレーションの最前線ではどんな研究が進んでいるんですか？
特に注目されている方向性を5つ紹介しよう：
- 物理インフォームドニューラルネットワーク（PINN）との融合：サロゲートモデルにPINNを使い、支配方程式の制約をサロゲート自体に組み込む。データ効率が飛躍的に向上する。
- マルチフィデリティキャリブレーション：粗メッシュの安価な計算と細メッシュの高精度な計算を階層的に組み合わせる。低フィデリティモデルで大域探索、高フィデリティモデルで精密化という戦略。
- デジタルツインとの統合：運用中のセンサデータをリアルタイムで取り込み、事後分布を逐次更新する。Bayesian filteringの技法を応用。
- 変分推論（VI）：MCMCの代替として事後分布を既知の分布族で近似する。収束が速いが、多峰性の見落としリスクあり。
- 確率的プログラミング言語：PyroやTensorFlow Probabilityなど、深層学習フレームワーク上でベイズモデルを構築する流れ。GPU加速でMCMCも高速化。
デジタルツインにベイズキャリブレーションを使うって、なんかSFっぽくてワクワクします！
実際、航空宇宙や原子力ではすでに実用段階に入っているよ。たとえばNASAのデジタルツインプログラムでは、飛行中のひずみゲージデータでFEMモデルの材料パラメータをリアルタイム更新している。日本でも原子力規制庁がV&Vの枠組みにベイズキャリブレーションを位置づけている。
トラブルシューティング
よくあるエラーと対策
ベイズキャリブレーションを走らせてみて「なんかおかしい」ってなったとき、まず何を確認すべきですか？
トラブルシューティングの定番チェックリストはこうだ：
| 症状 | 原因 | 対処法 |
|---|---|---|
| 採択率が極端に低い（<1%） | 提案分布のステップ幅が大きすぎる | 提案分布の共分散行列を縮小。AMやDRAMなど適応的手法に切り替え |
| 採択率が高すぎる（>95%） | ステップ幅が小さすぎて局所に停滞 | ステップ幅を拡大。目標採択率 23%（高次元）〜44%（1次元）を目安に |
| $\hat{R} \gg 1$ で収束しない | 多峰性、パラメータ間の強い相関、初期値が不適切 | チェーン数を増やす。異なる初期値で複数チェーン。TMCMC/SMCに切り替え |
| 事後分布が事前分布とほぼ同じ | データの情報量不足、パラメータの非識別性 | 感度分析で影響度を確認。実験データの追加。パラメータ数を削減 |
| 事後分布が非常に狭い（過信） | モデル不整合の無視。$\sigma$ が過小 | Kennedy-O'Haganの $\delta$ を導入。$\sigma$ を推定対象に含める |
| 対数尤度が $-\infty$ になる | パラメータがFEMの発散を引き起こす | 事前分布で物理的に安定な範囲に制約。FEM発散時の処理（尤度=0）を実装 |
| サロゲートの予測が不安定 | 訓練データが事後高確率領域をカバーしていない | 逐次的にFEM計算点を追加してサロゲートを再学習。Active learningの導入 |
まとめると、採択率・$\hat{R}$・トレースプロット・事後予測のデータ整合性、この4つを真っ先に確認すればいいですね？
その通り。あとは事後予測チェック（Posterior Predictive Check, PPC）も忘れずに。事後分布からパラメータをサンプリングしてモデルを走らせ、その予測分布が実験データを包含しているかを確認する。PPCで実験点が95%予測区間の外に多数出るなら、尤度関数かモデル自体を見直す必要がある。
ベイズキャリブレーション、奥が深いですけど体系的に理解できました！ まずはPyMCで簡単な例を試してみます。ありがとうございました！
いいね。PyMCの公式チュートリアルに線形回帰のベイズ推定があるから、まずそこから始めて、次にCAEモデルとの連携に進むのがおすすめだ。収束診断は絶対に省略しないこと！
関連トピック
なった
詳しく
報告