numpyでの円周率の計算（モンテカルロ法）

備忘録として簡単に書くことにする。

モンテカルロ法を用いた円周率の導出

import numpy as np


# 基準とする数(大きくするほど正確な円周率が求まる)
base_number = 100000
# 0〜1の一様乱数を２つ用意する
random.seed(0)
sample_x = np.random.uniform(0.0, 1.0, base_number)
random.seed(1)
sample_y = np.random.uniform(0.0, 1.0, base_number)

# 座標の原点からの距離の二乗の関数
def sample_function(x, y):
    return(x**2 + y**2)

# 半径1の円の中にいくつ点が存在しているか
count = 0
for x, y in zip(sample_x, sample_y):
    if sample_function(x, y) <= 1:
        count += 1
# 確認
print("count=", count)

# 以下のように考えられる
# (半径1の円を囲む一辺2の正方形の面積):(半径1の円の面積) = 4:π
# (半径1の円を囲む一辺2の正方形の面積):(半径1の円の面積) = (base_number):(count)
# つまり、xとyは一様乱数なので、全体の点の個数と円の内部に含まれる点の個数との比がおよそ4:πになるはず、ということ
# 以上より、
print("π=",  4*count/base_number)

出力結果

count= 78654
π= 3.14616

ある程度正確な値を得ることができた。なお、今回は円の内部に存在するかどうかを原点からの距離の二乗の値が１より小さいかどうかで判断した。もちろん、ルートで実際の距離にして判断しても良いがこのままで何も問題がないのでこのままにした。もしも実際の距離を出したいのであれば、

import math

math.hypot(x, y)

でxとyのそれぞれの二乗の和の正の平方根、つまり原点からの距離を計算できる。