pythonによるモンテカルロ法で円周率を近似する

本記事の内容は筆者の理解に基づいており、誤りが含まれる可能性があります。
内容に関してご指摘があれば、ぜひコメントでお知らせください。

モンテカルロ法とは

シミュレーションや数値計算を乱数を用いて行う手法の総称。
Wikipediaより引用

手順

１．0.0 ~ 1.0 の乱数x,yを２つ生成する。
２．座標(x,y)が半径１の円の内側か外側かを判定する。
３．πを求める。

プログラム

こちら

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

N = 100000
x = np.random.uniform(0.0,1.0,N)
y = np.random.uniform(0.0,1.0,N)

ins = []
out = []
count = 0
for i in range(N):
    if math.hypot(x[i],y[i])<1:
        ins.append((x[i],y[i]))
        count += 1
    else:
        out.append((x[i],y[i]))
ins_x,ins_y = zip(*ins)
out_x,out_y = zip(*out)
plt.figure(figsize=(8,7))
plt.plot(ins_x,ins_y,'o',label='Inside')
plt.plot(out_x,out_y,'o',label='Outside')
plt.legend()
plt.show()
print('円周率の近似値',4.0 * count / N)

解説

細かくやっていきます

とりあえず、必要な物をimport

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

乱数の生成

乱数は一様乱数です。
乱数を100000個生成させる。（多いほど精度は上がる）

一様乱数　すべての数が同じ確率で出現する乱数

N = 100000
x = np.random.uniform(0.0,1.0,N)
y = np.random.uniform(0.0,1.0,N)

半径1の円の内側か外側か判定

math.hypot(x,y)は　

\sqrt{x^2+y^2}

　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　　を求める。
math.hypot(x,y)<1なら内側　→　insリストに追加(inside)
そうでなければ外側 →　outリストに追加(outside)
count で円の内側にある点の数をカウント

ins = []
out = []
count = 0
for i in range(N):
    if math.hypot(x[i],y[i])<1:
        ins.append((x[i],y[i]))
        count += 1
    else:
        out.append((x[i],y[i]))

座標に分割

(x,y)はタプルなので、グラフにするために分割する

ins_x,ins_y = zip(*ins)
out_x,out_y = zip(*out)

グラフに描写

plt.figure：グラフサイズを設定
plt.plot(x座標,y座標) 、'o'は点で表示、labelはラベル（）
plt.legend()：ラベルを表示させる
plt.show()：グラフ表示

plt.figure(figsize=(8,7))
plt.plot(ins_x,ins_y,'o',label='Inside')
plt.plot(out_x,out_y,'o',label='Outside')
plt.legend()
plt.show()

グラフ

円周率の近似

グラフ(1.0 * 1.0)の正方形の面積　→　1
グラフの青（円の内側）の点の面積　→　π/4 (1.0 × 1.0 × π × 1/4 なので)
面積比
グラフの面積：扇形の面積 = 全体の点の個数：青の点の個数

　　　1　　　:　　\frac{π}{4}　　=　　　N　　　　:　　count

π = \frac{count}{N}　× 4

#小数点がつくので、わかりやすく4.0
print('円周率の近似値',4.0 * count / N)

そうすると、出力は約3.13ぐらいになると思います。
ぜひ、Nの値を変えてみましょう

最後に

これから、私のOutputとしていろいろ記事を投稿する予定です。

Thank you for reading!