Pythonで学ぶ高校数学の美しい物語（モンテカルロ法）

※このサイトの記事内では広告を掲載している場合があります。

今回は、Pythonを用いてモンテカルロ法によって円周率を求め、精度評価のシミュレーションをしていきます！楽しんで読んでください！

今回のコードの参考になると思われる記事をまとめておきます。

参考：【Python入門】初めてのプログラミング（ラムダ式 : 無名関数, lambda）
参考：【Python入門】初めてのプログラミング（スコープと関数）
参考：年金の最適受給月はいつか？Pythonでシミュレーション

プログラミングスクールに関しては下の記事で詳しく記述しています。

こちらで紹介しているスクールは、すべて無料期間がある優良なスクールのみで、特徴を明確にし、読者のニーズに絞って丁寧に解説しました。初めての一歩として、無料説明会に参加してみてください。

2021年9月25日【全スクール無料体験あり！】プログラミングを学びたい人におすすめのスクール5選！

モンテカルロ法とは

モンテカルロ法とは、乱数を用いて、シミュレーションや数値処理を行う手法です。

下の図を見てください。青い点は、乱数を用いて打った点だと考えてください。青い点を数えると、（左：右）＝（10：10）＝（1：1）になります。なので、面積比も１：１になります。

今回は、説明のため点の数は20個ですが、この点の数が10万以上あれば、より正確な値が得られ、誤差は0.1%以内になるといわれています。

円周率の求め方

今回は、円周率である、「π」の近似値を求めます。次のような条件の元考えていきます。

正方形の1辺の長さを2cmとすると、内接する円の半径は1cmとなる。
青い点の総数をNとする。
円の中にある青い点の数をnとする。

これらの条件と、面積比、円の面積を求める公式から、

2 × 2 ： 1 × 1 × π ＝ N ： n

よって、

π = 4n / N

となります。簡単ですね！

モンテカルロ法シミュレーション

ここからはPythonを用いてモンテカルロ法のシミュレーションを行っていきます。

プログラムの作成

単純に円周率を求めるだけなら以下のコードでできます。
何をしているのかが理解できれば大丈夫です。

円の原点を（0, 0）としています。


import random
import math

def Monte_Carlo():
    # 原点と任意の点までの距離を求める関数
    distance = lambda x, y : math.sqrt(x ** 2 + y ** 2)

    # 試行回数の設定
    N = int(input("試行回数"))

    # 円の中にあることを判定する変数
    R = 1

    # 円の中にある点を数える変数
    n = 0

    for _ in range(N):
        # -1 ~ 1までのランダムな浮動小数点（float型）
        x = random.uniform(-1, 1)
        y = random.uniform(-1, 1)

        r = distance(x, y)
        if r <= R:
            n += 1

    print(4 * n / N)

if __name__ == "__main__":
    Monte_Carlo()

"""
試行回数1000
3.216
"""

ただ、これだけだとつまらないのでコードを書き換えて視覚的にみていきましょう！

ラムダについては次の記事を参考にしてください。

参考：【Python入門】初めてのプログラミング（ラムダ式 : 無名関数, lambda）

10万回の視覚的試行

次は、どのように点が取られているのかを視覚的にわかるようにグラフ描写していきます。


import matplotlib.pyplot as plt
import numpy as np
import random
import math

def Monte_Carlo():
    distance = lambda x, y : math.sqrt(x ** 2 + y ** 2)

    # 試行回数
    N = 100000

    # グラフ描写用配列
    x_list = np.zeros(N)
    y_list = np.zeros(N)
    R = 1
    n = 0

    for i in range(N):
        x = random.uniform(-1, 1)
        x_list[i] = x

        y = random.uniform(-1, 1)
        y_list[i] = y

        r = distance(x, y)
        if r <= R:
            n += 1

    print("π :", 4 * n / N)
    print("内側の点 :", n, "個")
    print("外側の点 :", N - n, "個")

    draw_graph(x_list, y_list)

def draw_graph(x_list, y_list):
    c1 = plt.Circle((0, 0), radius=1, fc="None", ec="r", linewidth=2)
    ax = plt.gca()
    ax.add_patch(c1)
    plt.axis("scaled")
    plt.xlim(-1, 1)
    plt.ylim(-1, 1)
    plt.xlabel("X")
    plt.ylabel("Y")
    plt.scatter(x_list, y_list, marker=".", color = "blue", label = "POINT")
    plt.show()

if __name__ == "__main__":
    Monte_Carlo()

"""
π : 3.14472
内側の点 : 78618 個
外側の点 : 21382 個
"""