【Python/NumPy】モンテカルロ法の実装（円周率の計算）

Pythonでモンテカルロ法の実装し、円周率を計算する方法をソースコード付きでまとめました。

【モンテカルロ法】円周率の計算
【Matplotlib】グラフ化
【NumPy】モンテカルロ法の高速化

【モンテカルロ法】円周率の計算

モンテカルロ法とは、乱数を用いてシミュレーションや数値計算を行う手法の1つです。
モンテカルロ法で円周率を求める場合、以下の操作を行います。

手順	操作内容
①	正方形(1×1)内に点をランダムで生成します。
②	「生成した点」と「原点」の距離が1以下なら1ポイント、1より大きいなら0ポイントを与えます。
③	手順①②の操作をN回繰り返し、合計ポイントPを計算します。
④	合計ポイントPを4倍し、試行回数Nで割った値が円周率（近似値）となります。($\pi = \frac{4P}{N}$)
原理詳細	■【モンテカルロ法とは】円周率を計算するアルゴリスム

サンプルコード

サンプルプログラムのソースコードです。

# -*- coding: utf-8 -*-
import random
import time

def monte_method(N = 1000000):
    point = 0
    for i in range(N):
        x = random.random()
        y = random.random()
        if x*x+y*y < 1.0:
            point += 1

    # 円周率の近似解
    pi = 4.0 * point / N
    return pi


start_time = time.clock()

# 試行回数
N = 10000000

# モンテカルロ法
pi = monte_method(N)

# 実行時間を計算
end_time = time.clock()
run_time = end_time - start_time

# 実行結果
print("試行回数:", N) # 試行回数: 10000000
print("実行時間：", run_time) # 3.1853431999999997[sec]
print("円周率:", pi) # 円周率: 3.1418236

【Matplotlib】グラフ化

続いて、Matplotlibモジュールでグラフ化してみましょう。

# -*- coding: utf-8 -*-
import random
import matplotlib.pyplot as plt

def monte_method(N = 1000000):
    point = 0
    for i in range(N):
        # 乱数で点(x, y)をランダム生成
        x = random.random()
        y = random.random()

        # 原点からの距離が1未満（円内部）なら
        if x*x+y*y < 1.0:
            # ポイント加算
            point += 1

            # 赤点をプロット
            plt.plot(x, y,"ro")

        # 距離が1以上（円外部）なら
        else:
            # 青点をプロット
            plt.plot(x, y,"bo")

    # 円周率の近似解を返す
    return 4.0 * point / N

# 試行回数
N = 1000

# モンテカルロ法
pi = monte_method(N)

# 実行結果
print("試行回数:", N) # 試行回数: 100
print("円周率:", pi) # 円周率: 3.1418236

# 結果を表示
plt.grid(True)
plt.xlabel('X')
plt.ylabel('Y')
plt.show()

【NumPy】モンテカルロ法の高速化

先程の方法だと、for文を使っているため計算速度が遅くなります。
そこで、for文を使わずNumpy配列で一気に乱数の生成と個数のカウントを行います。

# -*- coding: utf-8 -*-
import random
import time
import numpy as np

def monte_method(N = 1000000):
    x = np.random.rand(N)
    y = np.random.rand(N)

    # 距離が1未満の点の個数をカウント
    point = np.sum(x*x+y*y < 1.0)

    # 円周率の近似解
    pi = 4.0 * point / N
    return pi


start_time = time.clock()

# 試行回数
N = 10000000

# モンテカルロ法
pi = monte_method(N)

# 実行時間を計算
end_time = time.clock()
run_time = end_time - start_time

# 実行結果
print("試行回数:", N) # 試行回数: 10000000
print("実行時間：", run_time) # 0.3774597[sec]
print("円周率:", pi) # 円周率: 3.1416828

実行時間が1/10になり、大幅に高速化できました。