計算物理学 No.8(モンテカルロ法 2)

統計力学、格子ゲージ理論など、さまざまな物理現象 (他に集団遺伝学、原子炉中の中性子拡散など) を表すのに広く使われるのがモンテカルロ法である。

1 マルコフ過程

x を物理系の状態を記述する 1 組の変数とする（各粒子の速度と運動量、イジングモデルでの全スピンの配置など）。数値シミュレーションにより一連の {x_{i}} が生成され、 {x_{i}} はある確率分布 P (x) に従う。

例えば、一様分布の乱数を使うと統計力学の物理量 f のカノニカル平均は

f_{n} = \frac{∑_{i=1}^{n} f(x_{i}) exp(−βH(x_{i}))}{∑_{i=1}^{n} exp(−βH(x_{i}))}

とした時

lim_{n→∞} f_{n} = ⟨f⟩

で与えられる。しかし、カノニカル平均では、平均エネルギー ⟨H⟩ 近くの状態が主に寄与するので、一様分布ではあまりにもサンプリングの効率が悪い。

そこで、あらかじめ確率分布を

P (x) ∝ exp(−βH(x))

ととることにすると、 f のカノニカル平均は

f_{n} = \frac{∑_{i=1}^{n} f(x_{i})}{∑_{i=1}^{n} １},  lim_{n→∞} f_{n} = ⟨f⟩

となり効率的である。

次の問題は、望む確率分布 P (x) に従う状態の系列 {x_{i}} を生成するルールを見つけることである。

確率分布 P (x) に従うシミュレーションの系列 {x_{i}} を導くため、状態 x_{i} がわかった時に、次の状態 x_{i+1} の確率分布を与えるルールを考察する。状態 x_{i} を与えた時に次の状態 x_{i+1} に移る遷移確率を W (x_{i}, x_{i+1}) と置く。この時、 W (x_{i}, x_{i+1}) が以下の条件を満たす時にマルコフ過程 (Markov process) と呼ぶ。 1. 規格化条件

∑_{y}W (x, y) = 1  (1)

2. エルゴード性

W (x, y) > 0  (2)

3. Limiting probability

∑_{x}P (x)W (x, y) = P (y)  (3)

[説明] 状態 x をベクトル空間とみなした時、遷移確率 W をその上の線形変換（行列）とみなすことが出来る。規格化条件とエルゴード性より、 W は（縮退のない）最大固有値 1 を持つ。 Limiting probability の条件より、最大固有値 1 に対応する固有ベクトルは P (x) であることがわかる。したがって、 p が十分大きい時に W^{p} を任意の状態に施すと、 P (x) の平衡分布に近付く。つまり、遷移確率 W にしたがって x_{i} を生成すると十分長い回数の後には望んだ確率分布 P (x) に近付く。

統計力学などでモンテカルロシミュレーションする場合は、遷移確率 W (x, y) をもっと単純な過程に分解し、

W (x, y) = ∑_{{z_{i}}}W_{1}(x, z_{1})W_{2}(z_{1}, z_{2}) ⋯ W_{k}(z_{k+1}, y)  (4)

とする (例えば状態 z_{j+1} は状態 z_{j} と局所的な違いしか無い)。この時、マルコフ過程の Limiting probability の条件は

P (x)W_{j}(x, y) = P (y)W_{j}(y, x)  (5)

という詳細釣合 (detailed balance) の条件に置き換えることが出来る。

[問題] (3) の左辺に (4)、 (5) を代入して (3) の右辺となることを確かめよ。

2 応用１：古典統計力学

2.1 熱浴法

統計力学で平衡状態を出すには

W_{j}(x, y) = P (y) ∝ exp(−βH(y))  (6)

としてみることが考えられる。

例として、イジングスピンの場合を考えよう。特定の位置にある１つのスピンのみ考え、残りを固定することで、次にその位置のスピンが上向きになるか、下向きになるかの確率を上の関係式に従い決めることが出来る。

この時、エネルギーの差を局所的に求めるだけで済むので（2 次元の平方格子では 1 つのスピンとその周辺スピン 4 つのみ）、比較的簡単な計算で済む。つまり、 j サイトでの局所場を

h_{j} = ∑_{k(j)} σ_{k}  (7)

(k(j) は j の隣接格子) としたとき、 j サイトのスピンが次に +1 となる確率を

p_{j} = \frac{exp(βh_{j})}{exp(βh_{j}) + exp(−βh_{j})} = (1 + exp(−2βh_{j}))^{−1}  (8)

−1 となる確率を

\frac{exp(−βh_{j})}{exp(βh_{j}) + exp(−βh_{j})} = 1 − p_{j}  (9)

とすれば良い。

2.2 メトロポリス法

自由度がイジングモデルのように離散的な場合は熱浴法で十分であるが、分子の運動の位置および速度のように連続自由度の場合は、局所的にも次のプロセスに無限の可能性が出るので、熱浴法では扱い難い。

この種の問題に対してはメトロポリス法が有効である。メトロポリス法の手続きは以下のとおりである。

1. 特定の箇所で、ランダムに試行状態 y を生成する。

2. H(y) ≤ H(x) ならば次の状態として y をとる。

3. H(y) > H(x) ならば、 [0, 1] の間の一様乱数 λ を作り、 λ < exp(−β(H(y) − H(x))) ならば次の状態として y をとる。それ以外なら以前の状態 x をそのまま維持する。

[問題]

メトロポリス法で長時間の確率分布がカノニカル分布になることを確かめよ。

参考文献

[1] Chapter 8 in “Statistical Field Theory” by Itzykson and Drouffe (Cambridge Univ. Press, 1989)