Pythonでモンテカルロ法を実装してみる

はじめに

どうも、初めまして。この度、Qiitaで書き始めたものです。どうぞよろしくお願いします。さて、最初の記事は何番煎じかはわかりませんが、モンテカルロ法をPythonで実装する方法を記事にすることにしました。

モンテカルロ法とは

そもそも、モンテカルロ法とはどんな方法か。ご存じの方も多いとは思いますが、今一度説明します。数値解析におけるモンテカルロ法は、確率を近似的に求める手法によく使われています。また、二次元上における閉曲面に囲まれた面積を求めることもできます。今回は、モンテカルロ法で円周率を求めてみます。

モンテカルロ法による円周率の求め方

まず、下の図のようなものがあるとします。

青線は半径が１の円の円周を表しています。ここに、ランダムで点を置いていき、置いた数（$n$）と円周内に置かれた点の数($p$)を数えます。そして、モンテカルロ法より以下の式が成り立ちます。

$$\frac{r^2\pi}{4S}=\frac{p}{n}$$

式内での$S$は正方形の面積です。$p$と$n$の比は面積の比と同じことから上の式が成り立つわけです。上の式から$\pi$について式変形すると、以下のようになります。

$$\pi=\frac{4Sp}{nr^2}$$

もちろん確率が絶対ではないため、あくまで近似値を求めることとなります。

Pythonで実装

ではパソコンで計算してみましょう。タイトルにもあるように、Pythonで実装してみました。以下はそのソースコードです。

pi_by_monte.py

import random

#変数の設定
in_pi = 0

#実行する回数
n = 100000

for i in range(n):
    x = random.uniform(0, 1)
    y = random.uniform(0, 1)

    z = x**2 + y**2

    if z <= 1:
        in_pi += 1

pi_by_monte = in_pi / n * 4

print(pi_by_monte)

意外と短いコードで実装出来ましたね。ソースコード内では実行回数を100000回としていますが、実際はもっと実行する必要があります。

結果

では、早速実行された結果を見てみましょう。

#結果はその時によって変わります
3.14885

とりあえず、3.14までは導出することに成功しました。しかし、それ以降は合っていません。上記したように確率が絶対ではないため、あくまでも近似値となるわけです。実行回数を1億回ほどに増やして実行してみました。

#結果はその時によって変わります
3.14180064

3.141まで合うようになりましたね。理論的にはこのまま実行回数を増やすともっと円周率に近くはなりますが、時間が非常にかかります。100000回ほどでしたらすぐに出ますが、1億回行うと3分以上もかかりました。このまま10億、100億と増やしていくと時間がどんどん増えていくのは目に見えています。もっと高速化する方法はあるとは思いますが、このプログラムではこれが限界ですね。ちなみに、10000回実行してmatplotlibを使用し、実際に点をプロットしたものを出してみました。（下の図）

円周よりも内側を赤色、外側の点を青でプロットしてみました。黒色の線が円周です。これがすべて埋まるくらいが理想ではありますが、メモリも相当食いますし処理するまでの時間がとんでもなくなるでしょうね。

まとめ

今回、Pythonを使いモンテカルロ法で円周率を求めてみました。簡単に求めることは出来ますが、精度や処理速度にやや難点ありと言ったとこではないでしょうか。しかし、実際にやって「モンテカルロ法ってこういうものなのか～」「モンテカルロスゲー」となるでしょう。ぜひ一度やってみてください。ご意見、質問などありましたらぜひコメントまでよろしくお願いします。