Matplotlib で円周率を求める

平成最後の日ですが、みなさん如何お過ごしでしょうか。私は家族旅行からもどり、明日5/1から心機一転、トライアスロンの練習と深層学習の勉強に集中したいと思ってる次第です。
　さて、前々から続けてた Matplotlib の学習で参考になったのが、乱数を用いて  1 x 1 の正方形の中にプロットして、円周率を求めるという方法。モンテカルロ法というらしいのですが、ロットする数をNとした場合、原点からの距離が1以内なら1ポイント、1以上なら0ポイントとして加算して総合ポイントをXとした場合、円周率は 4X / N で求められる。 

モンテカルロ法 - Wikipedia

シミュレーションや数値計算を乱数を用いて行う手法の総称。元々は、中性子が物質中を動き回る様子を探るためにスタニスワフ・ウラムが考案しジョン・フォン・ノイマンにより命名された手法。カジノで有名な国家モナコ公国の4つの地区（カルティ）の1つであるモンテカルロから名付けられた。ランダム法とも呼ばれる。

上記の計算を Matplotlib と numpyを使って計算して、実際に円周率と実行時間を求めてみよう。

import matplotlib.pyplot as plt
import numpy as np
import math
import time
%matplotlib inline

np.random.seed(100)
X = 0
N = 3000

fig_x = np.arange(0, 1, 0.001)
fig_y = np.sqrt(1- fig_x * fig_x)
plt.figure(figsize=(10,10))
plt.plot(fig_x, fig_y)

start_time = time.clock()

for i in range(0, N):
    #numpy.random.rand() で 0〜1 の一様乱数を生成する。
    score_x = np.random.rand()
    score_y = np.random.rand()
    if score_x * score_x + score_y * score_y < 1:
        plt.scatter(score_x, score_y, marker='o', color='b')
        X = X+1
    else:
        plt.scatter(score_x, score_y, marker='o', color='r')

pi = 4*float(X)/float(N)
end_time = time.clock()
time = end_time - start_time

print("円周率: %.10f"% pi)
print("実行時間:%f" % (time))

plt.grid(True)
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

実行結果

円周率 3.1506666667 と、まあまあ近似値が出た。乱数を用いて計算するという発想、この場合では地面に石を投げる方法でも求められるのだけど、数学者ってすごいですね。。そしてプログラミングで簡単に表現、計算できることも驚きです。