[Python] モンテカルロ法で円周率を求める

モンテカルロ法を用いて，円周率$\pi$を求めるプログラムのメモです．

例として，直径2cm（当然，半径1vm）の円の面積を求めることを考えます．
直径2cmの円は，1辺が2cmの正方形の中にぴったり入るので，円周率$\pi$を使えば面積は

$2 \times 2 \times \pi = 4\pi$

と求められます．

円周率$\pi$ がわからないときにどうすれば円の面積を求められるかと考えると，ある程度の点をランダムに打って，円の中に入った点と外れた点を数えて，その比率を求めてやれば，大体の円の面積が求められます（別に円でなくても構わない）．ランダムということを使って，何かを推定する手法のことをモンテカルロ法と言います．

以下のコードを実行してみます．

import random

import math

def monte_carlo(count):

    ok = 0

    for _ in range(count):

        x = random.random()

        y = random.random()

        if (x * x) + (y * y) <= 1:

            ok += 1

    return 4 * ok / count

print('For    100 {0}' .format(monte_carlo(100)))

print('For   1000 {0}' .format(monte_carlo(1000)))

print('For  10000 {0}' .format(monte_carlo(10000)))

print('For 100000 {0}' .format(monte_carlo(100000)))

（このコードは，人工知能はどのようにして 「名人」を超えたのか?に記載されていたものです）

monte_carloという関数の中でx座標とy座標をランダムに生成させてそれが半径1[cm]の円の中に入った場合の点の数を数えます．countは，試行回数で，何点のランダムの点を生成させるかを表す変数です．

上記のコードを実行してみると，試行回数（ランダムな点を生成させる回数）を増やすにつれて円周率$\pi = 3.141592...$に近づいていくことを確かめることができます．

$ python3 monte_carlo.py

For    100 3.2

For   1000 3.172

For  10000 3.142

For 100000 3.14368

当然，REPLで上記のプログラムを実行しても同じような（というのは乱数 - random - を用いているので全く同じにはならない）結果が得られるはずです．

$ python3

Python 3.6.7 (v3.6.7:6ec5cf24b7, Oct 20 2018, 03:02:14) 

[GCC 4.2.1 Compatible Apple LLVM 6.0 (clang-600.0.57)] on darwin

Type "help", "copyright", "credits" or "license" for more information.

>>> import random

>>> import math

>>> 

>>> 

>>> def monte_carlo(count):

...     ok = 0

...     for _ in range(count):

...             x = random.random()

...             y = random.random()

...             if (x * x) + (y * y) <= 1:

...                     ok += 1

...     return 4 * ok / count

... 

>>> print('For    100 {0}' .format(monte_carlo(100)))

For    100 3.12

>>> print('For   1000 {0}' .format(monte_carlo(1000)))

For   1000 3.124

>>> print('For  10000 {0}' .format(monte_carlo(10000)))

For  10000 3.1204

>>> print('For 100000 {0}' .format(monte_carlo(100000)))

For 100000 3.14728

>>> 

上記のコードをJupyter Notebookで実行した結果はGitHubで見ることができます．