単に円周率を使いたいなら
import math
print(math.py)
で良いんだけども，3.14がどうやって求まっているのか自分でやってみたい．
というわけで確率的モデルで円周率を求めていきたいと思います．
第一象限の適当なところに点を打っていく
まずx座標とy座標それぞれで，0から1のランダムな数値を定めて点を打ちます
import random
x, y = random.random(), random.random()
x^2+y^2が1以内かどうか確かめる
このxとyをそれぞれ2乗してお互いに足し合わせた数が1以内であれば，点は半径1の円の内側にあることになります．
number = x ** 2 + y ** 2 #これが1以内であれば円の内側
何回も繰り返して近似値を見る
これを何回も繰り返して，繰り返した数を分母，1以内だった回数×4(全ての象限で同じ結果になると仮定する)を分子とすると，円周率が算出されます．試しに1000回くらい繰り返してみましょう．
import random
incount = 0
def GenerateRandom():
x, y = random.random(), random.random()
number = x ** 2 + y ** 2
return number
iteration = 1000
for ite in range(iteration):
check = GenerateRandom()
if (check < 1):
incount += 1
quadrant = 4 #全象限
print(incount * quadrant / iteration)
3.2
うーん，微妙．数が少なかったようなので1000万回で試してみましょう．
# 5秒くらいかかります
3.141382
良きかな．
可視化してみる
数字が出ただけだとあんまり感動しなかったのでmatplotlibで可視化してみましょう．
まずは1000回バージョン
import random
import matplotlib.pyplot as plt
def GenerateRandom():
x, y = random.random(), random.random()
number = x ** 2 + y ** 2
return [number, x, y]
iteration = 1000
for ite in range(iteration):
check = GenerateRandom()
if (check[0] < 1):
plt.scatter(check[1], check[2], c = 'red', s = 10)
else:
plt.scatter(check[1], check[2], c='blue', s = 10)
plt.title('Monte Carlo Method')
plt.xlabel("x")
plt.ylabel("y")
plt.show()
10秒くらい待つと…
イマイチ分かりにくいので1000万回バージョンをみてみましょう．
import random
import matplotlib.pyplot as plt
def GenerateRandom():
x, y = random.random(), random.random()
number = x ** 2 + y ** 2
return [number, x, y]
iteration = 10000000
for ite in range(iteration):
check = GenerateRandom()
if (check[0] < 1):
plt.scatter(check[1], check[2], c = 'red', s = 10)
else:
plt.scatter(check[1], check[2], c='blue', s = 10)
plt.title('Monte Carlo Method')
plt.xlabel("x")
plt.ylabel("y")
plt.show()
これを実行したら一旦パソコンから離れて散歩でもしましょう．帰ってきた頃には処理が終わっているはずです．
グラフを正方形にすればより分かりやすいですが，綺麗に弧を描いて分布しているのが分かります．
これを計算したからといって何か生まれるわけでもありませんが，興味を持った方は試していただけると嬉しいです．
