モンテカルロ法の数値解法

モンテカルロ法は、確率分布に従うサンプルを作り、観測量の期待値を標本平均として推定する数値解法である。相互作用が強く、多自由度で解析解が得にくい系でも、統計力学の平均量や分布を直接評価できる枠組みである。

参考ドキュメント

• N. Metropolis et al., Equation of State Calculations by Fast Computing Machines, J. Chem. Phys. 21, 1087 (1953) https://bayes.wustl.edu/Manual/EquationOfState.pdf
• A. D. Sokal, Monte Carlo Methods in Statistical Mechanics: Foundations and New Algorithms (lecture notes / chapter) https://www.physics.ntua.gr/ph47/doc/sokal.ps
• 伊庭 幸人, モンテカルロ法入門（講義資料, PDF） https://www.ism.ac.jp/~iba/kougi_2006_ism/c20061.pdf

1. 定式化：積分を平均に置き換える

求めたい量が

$$I=\int f(x)p(x)dx$$

の形で書けるとき、$x_i\sim p(x)$ を $N$ 個生成すれば

$$\hat{I}=\frac{1}{N}\sum _{i=1}^{N}f(x_i)$$

で推定できる。誤差は概ね $O(N^{-1/2})$ で減るが、分散が大きい場合や自己相関が強い場合は収束が遅い。

統計力学では、状態 $x$（スピン配列、原子配置など）に対して

$$p(x)=\frac{1}{Z}\exp (-\beta E(x)),Z=\sum _x\exp (-\beta E(x))$$

が典型であり、熱平均は

$$⟨O⟩=\sum _{x}O(x)p(x)$$

として求める。

2. 重要度サンプリング

一様に $x$ を引くのではなく、寄与が大きい領域を重点的にサンプルして分散を下げるのが重要度サンプリングである。目的分布 $p(x)$ から直接引けない場合、後述の MCMC により「定常分布が $p(x)$ である連鎖」を作る。

3. MCMC（マルコフ連鎖モンテカルロ）で分布を生成する

3.1 要件

遷移確率 $T(x\to x^{\prime })$ を設計し、定常分布が $p(x)$ になることを保証する。十分条件として詳細釣り合い

$$p(x)T(x\to x^{\prime })=p(x^{\prime })T(x^{\prime }\to x)$$

を満たすように作ることが多い。さらに、状態空間を偏りなく到達できること（エルゴード性）が必要である。

3.2 Metropolis 法

提案分布 $q(x^{\prime }\mid x)$ で候補 $x^{\prime }$ を作り、受理確率

$$A(x\to x^{\prime })=min(1,\exp [-\beta (E(x^{\prime })-E(x))])$$

で受理する（提案が対称の場合）。提案が非対称なら Metropolis–Hastings として

$$A(x\to x^{\prime })=min(1,\frac{p(x^{\prime })q(x\mid x^{\prime })}{p(x)q(x^{\prime }\mid x)})$$

とする。

3.3 Gibbs サンプリング

部分変数 $x_i$ の条件付き分布 $p(x_i\mid x_{\mathrm{\neg }i})$ から直接引ける場合、受理棄却なしに更新できる。格子模型の局所更新やベイズ推定でよく現れる。

4. ムーブ（更新候補）の設計

同じ Metropolis 更新でも、提案（ムーブ）の作り方が効率を支配する。

代表例

• スピン系：単一スピン反転、クラスター反転
• 原子配置：局所変位、原子交換（swap）、欠陥の移動
• 合金・サイト占有：元素ラベルの入れ替え（半大正準、組成揺らぎの導入）
• 格子定数・体積：$NPT$ で体積（格子）ムーブを追加

受理率は高すぎても低すぎても効率が落ちるため、提案幅（変位量、回転角など）を調整して適度な受理率と自己相関の短さを両立させるのが定石である。

5. 収束判定と誤差評価

5.1 burn-in（熱化）と平衡性チェック

初期状態依存を避けるため、十分な burn-in を捨てる。平衡到達の目安として

• エネルギーや秩序変数の時系列が定常化している
• 初期条件を変えても同じ平均に落ちる
• 温度・圧力・化学ポテンシャルを変えたときに整合する を確認する。

5.2 自己相関と有効サンプル数

MCMC サンプルは相関を持つ。積分自己相関時間 ${\tau }_{\mathrm{int}}$ を導入すると、平均の不確かさは

$$\mathrm{SE}(\overline{O})\propto \sqrt{\frac{2{\tau }_{\mathrm{int}}}{N}}{\sigma }_O$$

の形で増大する。したがって、見かけの $N$ ではなく有効サンプル数

$$N_{\mathrm{eff}}\approx \frac{N}{2{\tau }_{\mathrm{int}}}$$

を意識する必要がある。

5.3 ブロッキング・ジャックナイフ・ブートストラップ

相関データの誤差を評価するには、系列をブロック平均して分散の飽和を見る方法が有効である。派生量（比、非線形変換、フィッティング由来の量）にはジャックナイフが扱いやすい。

6. 平均とゆらぎから熱力学量を得る

ボルツマン分布のサンプル列が得られれば、平均だけでなくゆらぎから応答関数を計算できる。

例

• 比熱：

$$C_V=\frac{⟨E^2⟩-⟨E{⟩}^2}{k_B{T}^2}$$

• 磁化感受率：

$$\chi =\frac{⟨M^2⟩-⟨M{⟩}^2}{k_{B}T}$$

• Binder 比（相転移の有限サイズ解析で用いることが多い）：

$$U_4=1-\frac{⟨M^4⟩}{3⟨M^2{⟩}^2}$$

7. 数値解法としての加速アルゴリズム

7.1 クラスター更新：臨界近傍の停滞を抑える

臨界点近傍では相関長が伸び、局所更新は極端に遅くなる（critical slowing down）。クラスター法は相関した自由度をまとめて更新し、自己相関を大幅に短縮する。

7.2 レプリカ交換

多峰性のエネルギー地形（スピングラス、欠陥配置、相共存など）では、複数温度のレプリカを同時に走らせ、配置を交換することで低温の停滞を緩和する。交換受理率は通常

$$A=min(1,\exp [({\beta }_i-{\beta }_j)(E_j-E_i)])$$

で与えられる。

7.3 拡張アンサンブルと再重み付け

一つの温度で集めたヒストグラムから近傍の温度の物理量を推定する再重み付け、あるいは密度状態 $g(E)$ を推定して温度を跨いだ熱力学量を得る方法（マルチカノニカル、Wang–Landau など）がある。これらは自由エネルギー障壁の横断を助ける。

8. kinetic Monte Carlo（kMC）：時間発展を数値的に積分する

平衡分布のサンプリングではなく、遷移率 $\{k_{\mu }\}$ に基づいて離散事象を発生させ、連続時間を進める方法が kinetic Monte Carlo である。次に起こる事象を確率 $k_{\mu }/\sum _{\nu }{k}_{\nu }$ で選び、時間刻みを

$$\mathrm{\Delta }t=-\frac{\ln u}{\sum _{\nu }{k}_{\nu }}(u\in (0,1))$$

で更新する。拡散、欠陥反応、析出・成長、表面吸着などの長時間スケールに適する。

9. 数値実装上のポイント

• 乱数生成器：品質の低い乱数は系統誤差を生むため、検証された生成器を用いる
• 再現性：seed、サンプル数、捨てステップ、ブロック化条件を記録する
• 並列：複数鎖（独立レプリカ）や温度レプリカ交換が自然に並列化しやすい。乱数ストリームの独立性（分割方法）に注意する

まとめ

• モンテカルロ法の数値解法は、(1) 目的分布の定義、(2) 更新則（詳細釣り合い＋エルゴード性）の設計、(3) ムーブの工夫による自己相関短縮、(4) 収束判定と相関を考慮した誤差評価、の組で成立する。
• 臨界近傍や多峰性地形では、クラスター法・レプリカ交換・拡張アンサンブルなどの加速法が本質になり、kMC は遷移率モデルに基づく時間発展を与える。