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

統計力学など、さまざまな物理現象 (他に集団遺伝学、原子炉中の中性子拡散など) を表すのに広く使われるのがモンテカルロ法である。高次元の積分など他に有効な数値計算法が無い場合にも有効である。

1 疑似乱数

まず、モンテカルロ法を使う前に乱数を発生させる数値的方法を考察する。厳密にはプログラムで生成される数列は決定論的なものであるから疑似乱数と呼ぶべきである。

シミュレーションに用いる疑似乱数列は以下の性質を持つことが要求される。

1. 統計的性質がよい。

2. 計算効率が高い。

3. 周期が長い。

4. 再現性がある。

1.1 線形合同法

漸化式

I_{j+1} = mod(aI_{j} + c, M)  (1)

により発生される数列 I_{1}, I_{2}, I_{3}, …, (0 ≤ I_{j} < M) を疑似乱数列とするものを線形合同法 (1948 年頃提案された) と呼ぶ。ここで、 M は正の整数で法と呼び、 a, c は 0 ≤ a, c < M という整数でそれぞれ乗数、増分と呼ぶ。 I_{j} を M で割る x = I_{j}/M と、 [0,1) 区間の一様乱数 0 ≤ x < 1 を得ることが出来る。

線形合同法により生成される疑似乱数列はいつかは元に戻り、その周期は M 以下である。パラメータ a, c の選び方により I_{j} の性質は大きく異なる。例えば、 a = 0 では I_{j} = c、 a = 1, c = 0 では I_{j} = I_{1}、 a = 1, c = 1 では周期 M ではあるが、とうてい乱数列とは呼べない。

[課題 1] 線形合同法による疑似乱数列を生成するプログラムを作成してみよ。 1.1.1 パラメータの推奨値

パラメータ a, c の選び方には以下のようなものが推奨されている (文献 [1] より引用 )。「M が 2 の累乗なら mod(a, 8) を 5 または 3 とし（5 の方が安全）、定数項 c を奇数と

する。このとき、周期はちょうど M となり、その 1 周期分には 0 から M - 1 までの整数が 1 個ずつ現れる。

高速にするには c = 0 とし、初期値 I_{1} を奇数にする。ただし、周期は M / 4 になる。」

1.1.2 線形合同法の欠点

線形合同法は簡便かつ有効な疑似乱数発生法ではあるが、いくつか欠点がある。

1. 合同法乱数は、上位の桁はランダムだが、下位の桁はランダムでない。

2. 合同法乱数は、 1 個ずつ使えばランダムだが、いくつか組にして使えばランダムでない。

性質 1. から説明すると、線形合同法による疑似乱数列は下位 k bit だけを見ると周期は高々 2^{k} である。

[課題 2] [課題 1] で作ったプログラムに対し、上記の性質を確認してみよ。複数のパラメータに対し試してみること。

性質 2. に関してはたとえば 0 から 15 までの乱数で、 (x, y) の座標を生成することを考える。ここで周期 16 の線形合同法で乱数を生成してみる。すると、生成される (x, y) の座標は 16 / 2 = 8 通りにしかならない。 (x, y) の座標は 16^{2} = 256 個もあるのに、 8 個の座標しか生成できないのです。 (x, y, z) の場合は 16^{3} = 4096 個の可能な座標があるのに、 16 個の座標しか生成できません。

[課題 3] [課題 1] で作ったプログラムに対し、上記の性質を確認してみよ。複数のパラメータに対し試してみること。 gnulpot の機能を使うと、 2 次元プロット、 3 次元プロット（splot を使う) で可視化し易くなる。

1.2 M 系列

線形合同法は簡便かつ有効な疑似乱数発生法ではあるが、「1 個ずつ使えばランダムだが、いくつか組にして使えばランダムでない」などの欠点がある。

線形合同法のほかには、 M 系列乱数 (M-sequence random numbers) というアルゴリズムが有名である。詳しくは述べないが、参考文献として [2] をあげておく。

1.3 Mersenne Twister (MT)

この他、 1996-7 年に提案された、 Mersenne Twister (MT) と呼ばれる疑似乱数発生法もある。参考文献 [3] 参照。

2 モンテカルロ法

乱数を使った無作為抽出を利用して数学や物理の問題を解こうという手法。例：円周率 π の計算正方領域 0 ≤ x < 1, 0 ≤ y < 1 内の点を乱数で多数発生させ、 1/4 円内

y ≤ \sqrt{1 − x^{2}}

に入る点の割合を調べる。この割合は 1 : π/4 になる。

参考文献

[1] 「C 言語による最新アルゴリズム事典」奥村晴彦著、技術評論社 (1991))

[2] http://www.ysr.net.it-chiba.ac.jp/data/rand/node6.html http://www.edu.cs.kobe-u.ac.jp/Algorithm/pdf/Note03.pdf

[3] http://www.math.keio.ac.jp/ matumoto/mt.html http://www.emit.jp/mt/mt.html