モンテカルロ法で円周率を求めたい

牛乳を毎日1リットル飲んでます。ko_ya346です。
モンテカルロ法で円周率を求められるそうなので、実際にやってみます。

なぜ？

統計の勉強をしていたらMCMC（マルコフ連鎖モンテカルロ法）に出くわしました。
こいつを理解して実装までしたいのですが、そもそもモンテカルロ法ってなんだっけ？となったので
やってみることにしました。

どうやって円周率を求めるの

以下の手順で考えます。

1. 1辺2rの正方形の中にぴったり入る半径rの円を考えます
2. この正方形の中に、ランダムに点を打っていきます
3. たくさん点を打つと、点の数はその領域の面積に近づいていきます
4. 以下の式により円周率を求めます

実際にやってみる

点の数を変えて実際に計算してみました。
↓結果

点の数が10000個超えると、ほぼ3.14くらいになりますね。
何度かやってみましたが同様の結果となったので、まあ成功ですね。
こんなシンプルな方法で円周率を求められるんですね！
中々面白かったです。

以下は使用したコード

def get_random_point(n: int) -> float:
    """
    一様分布[-1.0, 1.0] からなる乱数を返す 
    [0.0, 1.0)の乱数でも同様に求まりますが、
    円のほうがわかりやすい気がしたのでこの設定にしてます。
    """

    return np.random.uniform(-1.0, 1.0)

def is_in_circle(x: float, y: float, r=1) -> bool:
    """
    x, yが半径rの中にあるか判定

    Args:
        x, y: [-1.0, 1.0]
        r: 半径
    """
    return x**2 + y**2 < r**2

# 半径
R = 1 
# 点の数
n_points = [10, 100, 1000, 10000]

pi_lst = []
in_x = []
in_y = []
out_x = []
out_y = []


for n_point in n_points:

    cnt_into_circle = 0

    for i in range(n_point):
        x = get_random_point(n_point)
        y = get_random_point(n_point)
        if is_in_circle(x, y, R):
            cnt_into_circle += 1
            in_x.append(x)
            in_y.append(y)
        else:
            out_x.append(x)
            out_y.append(y)

    pi = cnt_into_circle / n_point * 4
    pi_lst.append(pi)

    plt.figure(figsize=(5, 5))
    plt.title(f"num point: {n_point}, pi: {pi}")
    plt.scatter(in_x, in_y)
    plt.scatter(out_x, out_y)
    plt.show()

print(pi_lst)

参考

okazukun5.com