몬테카를로법에 의해 π를 추정한다.

방법

사각형 안에 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회 해본다"는 것을 했다. 이것은 사각형 안에 치는 점의 수와 π의 시뮬레이션 횟수라는 두 개의 변수를 갖는다.

이 변수를 변경하여

사각형 안에 치는 점의 수를 늘린다.

π의 시뮬레이션 횟수 늘리기 => π의 추정값의 평균은 참 값에 가깝고 분산은 작아진다.

라는 것이 실현된다.