Pythonとモンテカルロ法による円周率近似

Pythonで機械学習を勉強する前に、まずはPythonになれなくちゃな！というわけでいろいろ調べたらモンテカルロ法とかいうのが出てきたんで試しにやってみた。

モンテカルロ法とは

「乱数を使って何かしらの値を確率的に求めようとする方法」らしいです。これだけではよくわからんけど、まあとにかく何かやってみたかったのでいろいろ調べてたら、円周率を求めようとしているところが結構あったので、とりあえずやってみようと思った。

円周率を近似するアルゴリズム

こちらのサイトを参考にさせてもらいました。

モンテカルロ法と円周率の近似計算 | 高校数学の美しい物語

要はxy座標系で乱数によって作った点（-1<x<1, -1<y<1）をいっぱいプロットして、その点がある正方形に内接する円の中に入るかを判定して、その個数の比から円周率を求めようってことらしい（間違ってたらすみません…）。

1辺の長さが1の正方形の面積と半径が1/2の円の面積の比は 1:π/4 、プロットする乱数の合計（正方形の内側に入る点の数）と円の内側に入る点の数をそれぞれM個、N個とすれば、その比は M:N になる。

面積の比と乱数の個数の比は等しくなるはずだから、1:π/4 = M:N が成り立ち、方程式をπについてとけば、

π = 4N/M

というふうになるはず。これを実装する。

Pythonによる実装

とりあえずやってみた。

1. #coding: utf-8
2. import numpy as np
3. def get_pi(n):
4.     x_list = np.array(np.random.random(n))
5.     y_list = np.array(np.random.random(n))
6.     r_list = (x_list**2 + y_list**2) ** 0.5
7.     cir_point = 0
8.     #三平方の定理から、円の中に入っているかどうかを判断する。
9.     for i in range(n):
10.         if r_list[i] <= 1:
11.             cir_point += 1
12.     square_point = n
13.     pi = 4*cir_point / square_point
14.     print("円内の点の数: ", cir_point)
15.     print("円外の点の数: ", square_point)
16.     print("近似した円周率(点が{}個の時): ".format(n), pi)
18. if __name__ == "__main__":
19.     n = int(input())    #プロットする点の数を決める
20.     get_pi(n)

結果

n=10^2 ~ 10^7までの結果は以下のようになった。

円内の点の数: 83
円外の点の数: 100
近似した円周率: 3.32

円内の点の数: 782
円外の点の数: 1000
近似した円周率(点が1000個の時): 3.128

円内の点の数: 7922
円外の点の数: 10000
近似した円周率(点が10000個の時): 3.1688

円内の点の数: 78333
円外の点の数: 100000
近似した円周率(点が100000個の時): 3.13332

円内の点の数: 785338
円外の点の数: 1000000
近似した円周率(点が1000000個の時): 3.141352

円内の点の数: 7855042
円外の点の数: 10000000
近似した円周率(点が10000000個の時): 3.1420168

点の数が増えていくにつれて円周率っぽい値に近づいていくのがわかる。大体10^5~7くらいの乱数をプロットすれば小数第２位まではうまく近似できた。

課題というかただちゃんと調べないで知ってる方法でやっただけなのだが、for文を使っているので処理速度がめちゃくちゃ遅い（10^8個だと１分ほどかかる）。まったくNumpyを使いこなせてないので、これから精進していきたい。