【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

【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

【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/10になり、大幅に高速化できました。