どうしても池の面積を調べたい時のモンテカルロ法 ~計算物理学I (朝倉書店)を参考にpythonを使って~

モンテカルロ積分の説明から簡単なものの実装まで

モンテカルロ積分とその他の積分法の違い

　普通、積分というと微小面積や微小体積を足し合わせて計算するイメージがあるかと思いますが、モンテカルロ積分はそれらとは大きく概念が異なる積分の手法です。モンテカルロ積分を簡潔に言えば「乱数を使用する」積分と言えるでしょう。

「乱数を使用する」とは

　あなたは牧場の中の柵で囲われた領域の面積($S_{Fence}$)を求めたいと考えました。しかし、あなたの手元には物差しもスケールも巻き尺もありません。ただ、あなたは牧場全体の面積($S_{Ranch}$)を知っています。さてこのような状況で、どのようにして柵の中の面積を求めることができるでしょうか。

　答えは「牧場内に適当に(一様にランダムに)たくさんの石を投げる」です。

　ランダムに投げられた石は牧場内の柵の内外のいずれかに落ちますが、投げ方が一様でかつ投げた石の数が十分なほど多ければ、柵の中に落ちる数($N_{in}$)と柵の外側に落ちる数($N_{out}$)は、それぞれの領域の面積に比例すると考えられます。

　すると、柵の中の面積($S_{Fence}$)は以下のように表すことができるでしょう。

S_{Fence} = S_{Ranch} \times \frac{N_{in}}{N_{in}+N_{out}}

　
　このような面積や体積の求め方をモンテカルロ積分といいます。プログラムの中で「投石」は乱数生成によって実施することとなります。したがって、モンテカルロ積分は乱数を使用した積分方法と言えます。

モンテカルロ積分の利点

　多次元積分において、その他の積分よりも高速かつ誤差が少ない結果を得ることができます。今回は1次元のみですが、次回は簡単な10次元積分を実装したいと思います。

モンテカルロ積分の簡単なものの実装

　さてモンテカルロ法について簡単に学習したところで、pythonを用いて簡単なモンテカルロ積分のプログラミングを作っていきます。
　具体的には計算物理学I(朝倉書店)のP.101からP.102の課題を私なりに作成しました。課題の内容は以下です。

1．一辺が2の正方形で囲まれた円形の池(r=1)を考える.
2. 池の面積の厳密な値は$\iint dA = \pi$という既知の値である.
3. $-1 \le r_i \le 1$ の乱数列を発生する
4. $i = 1$から$N$について$(x_i,y_i) = (r_{2i-1},r_{2i})$とする。
5. $x_i^2 + y_i^2 <1$なら$N_{pond}=N_{pond}+1$, その他は$N_{box}=N_{box}+1$とする.
6. (5.76)を用いて面積を計算し, そこから$\pi$を求める.

(5.76)とは先ほど示した$S_{Fence}$を求める式と同様のものです。
さて、上記問題を解決するためのコードを作成してみましたが、かなりシンプルなものになりました。

montecarlo.py

import numpy as np
import random

#投げる石の個数(size)を与えることで池の面積を返す式を定義
def MonteCarlo_circle(size):
    Npond = Nbox = 0.
    list_r = np.zeros(2*size)

    #石の落ちた座標を記録するための配列を作成
    list_xin = list_yin = list_xout = list_yout = np.array([])

    #-1≦r_i≦1の範囲の乱数列を生成
    for i0 in range(2*size):
        ran_Num = random.uniform(-1, 1)
        list_r[i0] = ran_Num
    
    #石の落ちた場所と個数を数える。
    for i1 in range(0,size):
        length = list_r[2*i1]**2 + list_r[2*i1+1]**2
        if length < 1:
            Npond += 1
            list_xin = np.append(list_xin, list_r[2*i1])
            list_yin = np.append(list_yin, list_r[2*i1+1])
        else :
            Nbox += 1
            list_xout = np.append(list_xout, list_r[2*i1])
            list_yout = np.append(list_yout, list_r[2*i1+1])
            
    #それぞれの石の数によって面積を計算する。4は池を囲む正方形の面積
    Apond = 4*Npond/(Npond + Nbox)

    #石の落ちた場所と面積の計算結果を返す。
    return list_xin, list_yin, list_xout, list_yout, Apond

　使用したパッケージはnumpyとrandomです。
　numpyは必ず必要というわけではありませんが行列やベクトルの計算に便利なので、今後の応用のために導入しました。(最初はnumpyのnumpy.random.rand()で乱数列を作成しようと考えていたのですが、生成される乱数rが0≦r＜1の範囲であり、問題の条件を厳密には満たさなかったのでやめたという裏話もあります……。参考1)
　randomは乱数列を生成するために使用しました。random.uniform(a, b)がa <= b であればa≦N≦bであるようなの浮動小数点数Nを返してくれるためです。参考2)

　その他は問題の通りの挙動としていますが、今回は池に落ちた時とそれ以外の時の石の座標を記録してプロットしてみることとします。

実行結果

　投げる石の数を5000個としたときの実行結果をプロットすると以下のようになりました。(赤い線が池の外周で、池の中に落ちた石を水色、外に落ちた石を緑色にしました。)

　ちゃんと、池の中に落ちたものとそうでないものを分類できていることが分かりますね。また、これはたまたまですが池の面積も厳密解の$\pi$とかなり近い値になりました。(実際は新たに乱数列を生成するたびに結果は変わります。石の数をより増やすことや、同じ投石数での計算結果を平均化することなどで結果の変動は小さくできると考えられます。)

　これで安心して池の面積が分かりますね。よかったよかった。また、何か間違いや勘違い等がありましたら教えていただければ幸いです。

参考書籍
R.H.Landau・他 (2018)『実践Pythonライブラリー 計算物理学I -数値計算の基礎/HPC/フーリエウェーブレット解析』 (小柳義夫・他訳) 朝倉書店