モンテカルロ法で円周率の近似値を求めるアレ
モンテカルロ法で円周率の近似値を求めることができます。
有名な、内部に円が描かれている正方形に点をたくさん打つアレです。
式はこんな感じ。
\pi = \frac{ 円の中に入った数 }{打った点の合計} \times 4
なんでこれで円周率の近似値が求まるのか？
初めて理解した時はけっこう感動しました。
やってることは単純で簡単なので、解説してみます。
解説
ステップ1
まず半径1の円と、その円に外接するような正方形（1辺の長さが2）を描きます。
ステップ2
次に「この正方形の中に無作為に点を打ったとき、点が円の中に入る確率」を考えます。
この確率は2通りの方法で求めることができます。
【単純な確率計算で求める】
確率 = \frac{ 円の中に入った数 }{打った点の合計} ・・・①
【面積の比率で求める】
無作為に点を打つので、正方形に対する円の面積の比率がそのまま確立と同じ値になります。
正方形の面積は「縦 × 横」なので 2 × 2 で 4。
円の面積は「半径 × 半径 × π」なので 1 × 1 × π で π。
つまり確率は下記で求められる。
確率 = \frac{円の面積}{正方形の面積} = \frac{\pi}{4} ・・・②
ステップ3
最後にステップ2で求めた2つの確率式をもとに式変形をします。
①と②はともに同じ確立を表した式なので、「=」で繋げられます。
\frac{\pi}{4} = \frac{ 円の中に入った数 }{打った点の合計}
はい、これで両辺に「4」をかければ最初に紹介した式の完成です。
\pi = \frac{ 円の中に入った数 }{打った点の合計} \times 4
Pythonで実装
仕組みを理解したところで、実際に近似値が求まるかPythonで実装してみます。
円の中に入ったかどうかの判定には円の公式を利用します。
半径1の円の公式はこちら。
x^2 + y^2 = 1
つまりランダムに生成したx座標、y座標が以下の条件を満たせば、円の中に入ったと言えます。
x^2 + y^2 <= 1
では、点を100回打ってみます。
ソースコードはこちら。
import random
def get_pi(iteration):
in_square = 0
in_circle = 0
for i in range(iteration):
x = random.random()
y = random.random()
in_square += 1
if (x * x + y * y <= 1):
in_circle += 1
pi = (in_circle / in_square) * 4
return(pi)
print(get_pi(100))
出力結果：3.2
おお、円周率は「3.1415926535...」なのでまあまあ良い値が出ましたね。
このように、たくさんのランダムな値（乱数）を与えて、事象を確立的に解析することを
モンテカルロ法といいます。
試行回数を増やしてみる
先ほど点を打った回数は100回でしたが、もっと試行回数を増やしてみましょう。
| 試行回数 | 近似値 |
|---|---|
| 100回 | 3.2 |
| 1,000回 | 3.22 |
| 10,000回 | 3.1536 |
| 100,000回 | 3.13652 |
| 1,000,000回 | 3.148724 |
| 10,000,000回 | 3.1414596 |
| 100,000,000回 | 3.141593 |
試行回数を増やすほど精度が向上していますね。
せっかくなので可視化も
せっかくなので可視化もしてみました。
点を打つ回数が100回のものと10,000回のもので、それぞれ100個の近似値を求めて比較してみます。
点を100回打ったときよりも、1万回打って求めた近似値の方が「3.14...」に近い値が多いことが分かります。
可視化すると説得力と安心感が増しますね。
可視化大事。
ちなみに可視化はPythonのmatplotlibというライブラリを使用して行いました。
import matplotlib.pyplot as plt
x = []
for i in range(100):
x.append(get_pi(10000))
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.hist(x, bins=10, rwidth=0.8)
ax.set_xlim(2.7,3.6)
ax.set_ylim(0,30)
fig.show()
おわり。