モンテカルロ法 Python

monte.py

#! /usr/bin/python
# coding: utf-8
from random import *
from math   import *

# アークコサイン 定義 => y = cos**(-1)x |x| <=1, |y| <= pi / 2
# よくわかんなけど、
# アークコサインの角度(ラジアン）が-1のときは
# 円周率らしいので、比較用にそれなりの正確な近似値を出しておく
math_pi = acos( -1 )

seed()			# 乱数の種の初期化 今回はシステム時刻で乱数生成
whole = 2 ** 20	# 乱数発生回数 2の20乗 -> 1,048,576 だいたい１００万
k = 0.0			# 半径1未満の件数を格納する
for i in range(whole):
	x = random() # 値域[0.0, 1.0)のランダムな数を生成
	y = random() # 値域[0.0, 1.0)のランダムな数を生成
	if ( x**2 + y**2 < 1 ):
		k += 1
# 近似値をモンテカルロ法により出す
pi = 4 * k / whole
# 答え合わせ
print '乱数を使っているので、計算結果は毎度ちがうよー'
print '本当の答え               =' + str( math_pi )
print 'モンテカルロ法による答え =' + str( pi )
print '誤差                     =' + str( math_pi - pi )