モンテカルロ法における円周率の近似値の計算

monte_carlo.py

# -*- coding:utf-8 -*-
import numpy
import pylab
import math
import time 

X = 0 # 的に当たった回数
N = 10000 # 試行回数

# 四分円の境界の方程式[y=√1-x^2 (0<=x<=1)]を描画
circle_x = numpy.arange(0,1,0.001)
circle_y = numpy.sqrt(1- circle_x * circle_x)
pylab.plot(circle_x, circle_y) 

# N回の試行にかかる時間を計測
start_time = time.clock() 

# N回の試行を開始
for i in range(1, N):
    score_x = numpy.random.rand()
    score_y = numpy.random.rand()
    if score_x * score_x + score_y * score_y < 1:
        #的に入ったものは赤で表示
        pylab.plot(score_x, score_y,"ro")
        X = X + 1
    else:
        #的から外れたものは青で表示
        pylab.plot(score_x, score_y,"bo")

# piの近似値をここで計算する
pi = 4*float(X)/float(N)

# モンテカルロ法の実行時間を計算
end_time = time.clock() 
time = end_time - start_time

#　実行結果をUNIX上に表示する
print u"円周率の近似値:%f" % (pi)
print u"実行時間:%f" % (time)

# 結果を表示 
pylab.grid(True)
pylab.xlabel('X')
pylab.ylabel('Y')
pylab.show()