【Python実装】モンテカルロ法

本記事のアルゴリズムの記事はこちらです。

アルゴリズムが単純なので、ソースを貼り付けるだけにします。

###############################
#       モンテカルロ法
###############################
import numpy as np
import matplotlib.pyplot as plt

# 乱数を固定
np.random.seed(1)

# 円作成
x, y = [],[]

for d in np.linspace(-180, 180, 360):
    x.append(np.sin(np.deg2rad(d))) # sinの値をxとする
    y.append(np.cos(np.deg2rad(d))) # cosの値をyとする

# 点を発生
N = 100   # 発生させる点の個数

xmc = np.random.uniform(-1, 1, N)   # -1 から 1 まで N個乱数を発生させる
ymc = np.random.uniform(-1, 1, N)   # -1 から 1 まで N個乱数を発生させる
r = (xmc ** 2 + ymc ** 2) ** 0.5    # 原点からの距離を計算
accept = np.where(r <= 1, 1, 0)     # 円の中に点を1, 円の外にある点を0とする
accept_ratio = np.sum(accept) / N   # 円の中に点がある確率を求める
print("発生させた点の数: ", N)
print("数値的に計算 π: ", accept_ratio * 4)
print("真の値 π: ", np.pi)

# 円描画
plt.gca().set_aspect('equal', adjustable='box') # グラフを正方形にする
plt.plot(x,y)   # 円の線分をプロット

# 点描画
plt.scatter(xmc, ymc, color = "red", marker = ".")  # 点をプロット
plt.show()  # グラフを描画

実行結果は以下です。

おまけ

プログラムをちょっと変えて、を増やして数値的に求めたが真の値に近づくか調べてみました。
グラフの横軸はで、対数を取っていることに注意してください。

点の数		の値
10	1.0	3.2
50	1.69	2.64
100	2.0	2.88
500	2.69	3.168
1000	3.0	3.08
5000	3.69	3.1552
10000	4.0	3.1188
50000	4.69	3.12816
100000	5.0	3.1474
500000	5.69	3.14216
1000000	6.0	3.140636
5000000	6.69	3.141756
10000000	7.0	3.1407636

偉人の名言

どんなに悔いても過去は変わらない。どれほど心配したところで未来もどうなるものでもない。
いま、現在に最善を尽くすことである。
松下幸之助

動画