🚀
今更ながらモンテカルロ法の再勉強 ~円周率取得編~

2025/04/30に公開
そういえば、モンテカルロ法で円周率の取得ができるというのはずっと前からわかっていたもののやった試しがないなということで、今回は試してみようと思います。

 モンテカルロ法とはWikipediaによると、
モンテカルロ法（モンテカルロほう、（英: Monte Carlo method、MC）とはシミュレーションや数値計算を乱数を用いて行う手法の総称。元々は、中性子が物質中を動き回る様子を探るためにスタニスワフ・ウラムが考案しジョン・フォン・ノイマンにより命名された手法。カジノで有名な国家モナコ公国の4つの地区（カルティ）の1つであるモンテカルロから名付けられた。ランダム法とも呼ばれる。
ということです。とりあえず、乱数を用いて何かしらすることの総称とおもっておきましょう。それで間違いないはずです。
https://ja.wikipedia.org/wiki/モンテカルロ法

 モンテカルロ法でどうやって円周率を取得するのか？ではどうやってモンテカルロ法を利用すれば円周率を計算できるか説明していきます。
まずは、以下のように一辺の長さが2の正方形の中に半径1の円が内接している状態を考えます。
ここで、正方形の中にランダムな場所に点を打つことを考えます。この時、打たれた点が正方形の中に入っており、かつ円の中にも含まれている確率を計算してみます。点は必ず正方形の中含まれるので、母数は正方形の中に含まれている点の数（つまり全ての点の数）と等しくなり、これを踏まえるとこの確率は以下のような式で表せます。
円の中に点が入る確率 = \dfrac{円の中に点が打たれた数}{全ての点の数}また、この計算を図形の面積比で考えると以下のようにも計算できます。
円の中に点が入る確率 = \dfrac{円の面積}{正方形の面積} = \dfrac{1 \times 1 \times \pi}{2 \times 2} = \dfrac{\pi}{4}この二つの式を組み合わせると、
 \dfrac{円の中に点が打たれた数}{全ての点の数} = \dfrac{\pi}{4} となり、ここから\piは
 \pi = \dfrac{円の中に点が打たれた数 \times 4}{全ての点の数}と表せることがわかります。

 Pythonで検証してみるそれでは、Pythonを使ってこれが本当に近似値として計算できるか検証してみましょう。今回は点の数を増やした時に計算された結果についてグラフ化してみます。
X座標、Y座標を[-1, 1]の範囲で乱数により生成する
原点(0, 0）からの距離が1以下の点をカウントする
点を順番に100万回打ち、各タイミングでの計算結果をプロットする
全部で5回同じ条件で試行を行う
from math import pi
from tqdm import tqdm
import matplotlib.pyplot as plt
from random import uniform

TRIAL = 5
ITERATION = 1_000_000


def generate_coordinate() -> tuple[float, float]:
    return uniform(-1, 1), uniform(-1, 1)


def main():
    result = {}
    for trial in tqdm(range(TRIAL)):
        dot_in_circle = 0
        result[trial] = []

        for i in tqdm(range(1, ITERATION+1), leave=False):
            x, y = generate_coordinate()
            length_from_origin = (x ** 2 + y ** 2) ** 0.5

            if length_from_origin <= 1.0:
                dot_in_circle += 1

            result[trial].append(dot_in_circle * 4 / i)

        plt.plot(result[trial], label=f"trial {trial}")

    plt.xlabel("the number of dots")
    plt.ylabel("Calculated pi")
    plt.xscale("log")
    plt.hlines([pi], 0, ITERATION, linestyles="dashed")
    plt.grid()
    plt.legend()
    plt.show()


if __name__ == "__main__":
    main()
この実験結果は以下のようになります（横軸は値が大きいのでlogスケールとしています）。結果から、最終的には円周率におおよそ収束しているのではないかなと思います。なお、最終的な計算誤差は
{0: 0.0023983507257350927, 1: 0.016403985717814686, 2: 0.10056511962481182, 3: 0.08978419231309806, 4: 0.027735817665958368}
となっていたので、誤差は1%にも満たないということで、結構な精度なのではと思います。

 まとめ今回はモンテカルロ法を使って円周率を計算してみました。円周率の求め方としては結構有名な気がしますが、その他の方法もありますので、機会があれば試してみたいと思います。

 おまけみなさん、ラマヌジャンという天才数学者をご存知でしょうか？ここでは深くは触れませんが、名前を調べてもらうととてもたくさんの興味深い話が出てくるはずです。ここでは、彼が提案した計算式の一つを試してみたいと思います。その数式は以下です。
 \dfrac{1}{\pi} = \dfrac{\sqrt{8}}{99^2} \sum^{\infty}_{n=0} \dfrac{(4n)!}{\left( 4^nn!\right)^4} \dfrac{1103+26390n}{99^{4n}}言いたいことはわかります。私も同じ考えです。まずこれを計算するためのコードを書くのも嫌になるくらいですが、実装してみまして、intをfloatに変換する時にオーバーフローが起こったということでエラーになりまして、、、

NumPyやScipyを使ってやればいけると思いますが、実行しても階乗の計算もありますしパソコンに負担をかけたくなくて諦めてしまいました。もしやる気のある方は以下のコードをちょっと改修して試してみてください。
※機会があれば私の方で対応して記事更新またはコメントで共有します
from math import pi, factorial
from tqdm import tqdm
import matplotlib.pyplot as plt
from random import uniform

TRIAL = 5
ITERATION = 1_000_000



def calculate_error(calculated_pi) -> float:
    return abs(pi - calculated_pi) / pi * 100.0


def calculate_raj(i, prev):
    return (8**0.5) / (99**2) * (prev + factorial(4*i) / ((4 ** i) * factorial(i)) ** 4 * (1103 + 26390 * i) / (99 ** (4 * i)) )


def main():
    result = {}
    errors = {}

    for trial in tqdm(range(TRIAL)):
        dot_in_circle = 0
        result[trial] = []
        prev = 0


        for i in tqdm(range(1, ITERATION+1), leave=False):
            prev = calculate_raj(i, prev)
            result[trial].append(prev)

        plt.plot(result[trial], label=f"trial {trial}")
        errors[trial] = calculate_error(result[trial][-1])


    plt.xlabel("the number of iteration")
    plt.ylabel("Calculated pi")
    plt.xscale("log")
    plt.hlines([pi], 0, ITERATION, linestyles="dashed")
    plt.grid()
    plt.legend()
    plt.show()
    print(errors)


if __name__ == "__main__":
    main()
モンテカルロ法とは

モンテカルロ法でどうやって円周率を取得するのか？

Pythonで検証してみる

まとめ

おまけ

Discussion
