モンテカルロ法によってπを推定する。
方法
正方形の中に、x,y座標ともに一様な乱数を点としてプロットする。そのうち、円の中にある点の数を数え上げることで、円の中に存在する点の確率が計算できる。
この確率から面積を求めることができる。
具体的には、x,y座標ともに(-1,-1)と(1,1)となる正方形とそこに密接する円を想定する。
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
fig,ax = plt.subplots(figsize=(5, 5))
rect = patches.Rectangle((-1,-1),2,2,linewidth=1, fill=False, color='blue')
ax.add_patch(rect)
ax.add_artist(plt.Circle((0, 0), 1, fill=False, color='r'))
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(-2,2)
ax.set_ylim(-2,2)
この赤の円の中に入る点の数を数える。
1000個の点を打ってみる
df = pd.DataFrame(
(1 + 1) * np.random.rand(1000, 2) - 1,
columns=['x', 'y'])
fig, ax = plt.subplots(figsize=(5, 5))
ax.scatter(df.x, df.y, s=1, color='black')
rect = patches.Rectangle((-1,-1),2,2,linewidth=1, fill=False, color='blue')
ax.add_patch(rect)
circle = plt.Circle((0, 0), 1, fill=False, color='r')
ax.add_artist(circle)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_xlim(-2,2)
ax.set_ylim(-2,2)
この点群のうち幾つが円の中にあるかを計算する。
df["is_inside_circle"] = df.x**2 + df.y**2 <= 1
len(df[df.is_inside_circle]) # 798
数え方は
$$ x^{2} + y^{2} \leq 1 $$
となる点の数を数えている。
結果として、1000個のうち、798個が円の中に存在する点である。
正方形の中の79.8%が点が円を構成していると考えると、円の面積は
$$ 2 \times 2 \times 0.798 = 3.192 $$
円の公式では、π=3.14とすると
$$ 1 \times 1 \times 3.14 = 3.14 $$
より近い値が出ている。
1000個の点の配置の仕方は、numpyの乱数の生成の仕方に由来する。そのため、円の中にたまたま点が多く置かれるような偏りが起こる。
そこで、試行回数（シミュレーション回数）を増やすことで、計算されるπの値の確率分布を探ることができる。
それでは、これを1000回シミュレーションしてみる。
四角形の中に点を1000個打ち、πを計算する作業を1000回やってみる
results = [] # この変数のなかに推定されるπの値を入れる
for i in range(1000):
rand_num = 1000
xy = (1 + 1) * np.random.rand(rand_num, 2) - 1
df = pd.DataFrame(xy, columns=['x', 'y'])
df["is_inside_circle"] = df.x**2 + df.y**2 <= 1
ans = 4*len(df[df.is_inside_circle == True])/rand_num
results.append(ans)
# 描画のために整形
r_df = pd.DataFrame(results, columns=['area'])
plt.hist(r_df.area, bins=50)
print(f"mean: {np.mean(results)}") # 平均: 3.1412359999999997
print(f"variance: {np.var(results)}") # 分散: 0.0026137523040000036
1000回試行してみることで
「1000個の点を用いたモンテカルロ法では、πは平均3.1412、分散0.0026の範囲の値を取りうる」ことがわかった。
分散の平方根を取ることで、標準偏差がわかる。
np.sqrt(np.var(results)) # 0.051124869721105436
πの推定値の確率分布が正規分布に従っていると仮定すると以下のことがわかる。
値の68%は平均から±1σ、値の95%は平均から±2σ、99.7%は平均から±3σ以内に入る。（σ = 標準偏差 = ここでは0.051）
実際に計算してみよう。
平均から±1σの範囲にある点の個数
result_mean = np.mean(results)
sigma = np.sqrt(np.var(results))
# 平均から±1σの範囲にある点の個数
len([result for result in results if result >result_mean + (-1) * sigma and result < result_mean + 1 * sigma])
# 686
# 平均から±2σの範囲にある点の個数
len([result for result in results if result >result_mean + (-2) * sigma and result < result_mean + 2 * sigma])
# 958
# 平均から±3σの範囲にある点の個数
len([result for result in results if result >result_mean + (-3) * sigma and result < result_mean + 3 * sigma])
# 998
より、正規分布に近いような個数が得られることが分かる。
まとめ
この記事では、「四角形の中に点を1000個打ち、πを計算する作業を1000回やってみる」というのをおこなった。これは、四角形の中に打つ点の数と、πのシミュレーション回数という2つの変数を持つ。
この変数を変えることで
- 四角形の中に打つ点の数を増やす => 一試行あたりのπの値が真の値に近くなる。
- πのシミュレーション回数を増やす => πの推定値の平均は真の値に近くなり、分散は小さくなる。
ということが実現される。