5 モンテカルロ積分法

5.1 単純な場合(2次元の面積)

ここでは,モンテカルロ法(Monte Carlo method)呼ばれる乱数を使った一風変わった数値 積分の方法を示す.モンテカルロとは有名なモナコのカジノで,博打をするところである. この方法を少し学習すると,名前の由来が分かったような気になるはずである..

まずは,図5の様な円の面積を求めることにする.もちろん,この 円の面積は$ \pi$と分かっているが,これを乱数を用いて数値計算しようと言うのである. 円を囲む正方形内部には,ランダムな点がある.正方形内部にある点の個数を$ N_s$,円の 内部にある個数を$ N_c$とする.そして,正方形の面積を$ S_s$,円を$ S_c$とする.する と,これらには

$\displaystyle \frac{S_c}{S_s}\simeq\frac{N_c}{N_s}$ (28)

の関係がある.ランダムな点の数増加させると,これらは,ほとんど等式で結ばれるであ ろう.形状が単純な正方形の面積$ S_s$は分かっているので,円の面積$ S_c$

$\displaystyle S_c\simeq\frac{N_c}{N_s}S_s$ (29)

と計算できる.点の数を増やすと,それがランダムに分布する限り,右辺は円の面積に近 づく.正方形の面積と点の数を計算するだけで円の面積が分かるのである.このようにラ ンダムな点(乱数)を用いて,積分する方法をモンテカルロ積分と言う.

円の面積を求めても大したことはないが,例えば図6のよ うな複雑な形状の面積が計算できるとなるとぐっと御利益はある.先ほどの円の面積と同 じようにして,計算できるのである.計算したい部分内の点の数$ N_{in}$は,

$\displaystyle N_{in}$ $\displaystyle =($大円$\displaystyle \Omega_0$内部の点数$\displaystyle ) -($大円$\displaystyle \Omega_0$内部かつ小円$\displaystyle \Omega_1$内部の点数$\displaystyle )$    
  $\displaystyle \qquad\qquad -($大円$\displaystyle \Omega_0$内部かつ楕円$\displaystyle \Omega_2$内部の点数$\displaystyle ) -($大円$\displaystyle \Omega_0$内部かつ正方形$\displaystyle \Omega_3$内部の点数$\displaystyle )$ (30)

こうやって,内部の点の数を調べて,後は,式(30)のようにすれ ば面積を求めることができる.

このように複雑な領域の積分でも適用でき,応用範囲の広い方法である.積分の 精度を上げようとすると,サンプル数$ N$を増やすことになる.しかし,その精度は,た かだかサンプル数の平方根に比例するのみである.誤差を1/10にしたければ,サンプル数 を100倍にする必要がある.計算時間も100倍になる.これが,この方法の最大の弱点である.

精度がサンプル数の平方根に比例するのは,完全な乱数を仮定しているからである.それ を計算機で発生させることは不可能なので,その生成も問題となる.できるだけ良い乱数 を使わなくてはならないが,この問題は奥が深い.ここではこれ以上立ち入らないが,興 味のある者は調べてみると良い.

話は変わるが,ひとつ注意を与えておく.実際のプログラムでは,乱数を使って$ (x,y)$ を発生させて,それが内部にあるか否かを即座に判断する.配列に座標データを入れて, 最後に検査するとメモリーがいくらあっても足りなくなるからである.

図 5: モンテカルロ積分.正方形内にランダムな点をばらまく.ランダムな点の数 を増やせば,正方形内と円内の点の数の比が,正方形と円の面積の比になる.
\includegraphics[keepaspectratio, scale=1.0]{figure/MonteCarlo.eps}
図 6: この面積を解析的に計算するのは,大変面倒である.
\includegraphics[keepaspectratio, scale=1.0]{figure/MonteCarlo_complex.eps}

5.2 複雑な積分

もう少し複雑な積分が必要な場合を考えてみよう.ある関数 $ f(x_1,x_2,x_3,\cdot,x_M)$ の体積積分を考える.この体積積分は

$\displaystyle \int_\Omega f(x_1,x_2,x_3,\cdots,x_M)dx_1dx_2dx_3\cdots dx_M= V\langle f \rangle$ (31)

となる.ここで,$ V$はM次元体の領域$ \Omega$の体積, $ \langle f \rangle$はその内部 の関数の平均値である.これは3次元体の質量を考えると簡単である.その質量は,密度 の体積積分となる.これは体積に平均密度を乗じた値に等しい.当たり前の式である.こ のことから,式(32)の積分の値が欲しければ,体積と平均値が分か ればよいことになる.

積分を計算するために,まずは体積である.先ほどは,二次元問題であるが,はかなり複 雑な形状の面積をモンテカルロ積分で計算する方法を示した.2次元と同じ考えで,3次元, 4次元,$ \cdots$,M次元の体積も比較的,容易に求めることができる.分かり切った体積 に,領域$ \Omega$を包み込み,その内部にランダムに配置されたサンプル点の数を数えれ ば良いのである.体積の計算は簡単である.

残りは,体積内部の平均 $ \langle f \rangle$である.これも簡単で,領域$ \Omega$内部 にあるサンプル点の平均より求めることができる.即ち,

$\displaystyle \langle f \rangle\simeq\frac{1}{N_{in}}\sum_i f(\boldsymbol{r}_i)$   領域$ \Omega$の内部のみ (32)

となる.ここで, $ \boldsymbol{r}_i$$ i$番目のサンプル点の $ (x_1,x_2,x_3,\cdot,x_N)$座標で ある.また,$ N_{in}$は領域内部のサンプル数である.この計算も簡単で,コンピューター は難なく,平均値に近似値を求めることができる.

以上より,モンテカルロ法を用いると,体積$ V$と平均値 $ \langle f \rangle$の近似値が 計算できることが分かる.従って,式(32)の近似値を求めることが できる.

数値計算で有名な本「NUMERICAL RECIPES in C」 [2]には,図 7の様な形状の重心を求める問題が出ている.このように複雑な積 分が必要なときにモンテカルロ法は威力を発揮する.

図 7: 「NUMERICAL RECIPES in C」の例題.一部を切り取られたトーラスの重心を 求める.
\includegraphics[keepaspectratio, scale=0.5]{figure/NRC_example.eps}

ホームページ: Yamamoto's laboratory
著者: 山本昌志
Yamamoto Masashi
平成19年1月19日


no counter