モンテカルロ法を使ってみる（円周率の計算）

お遊びです。

乱数を使って数値計算やシミュレーションをおこなう
モンテカルロ法を使って円周率を求めてみた、という記事です。
シンプルでわかりやすく実装もしやすいので、
色んな人が解説したりブログ記事書いたりしている印象があります。

環境はGoogle Colaboratory(python3, jupyter notebook)です。

目次

1. モンテカルロ法
2. 実装＆計算結果
3. まとめ

1. モンテカルロ法

名前のもととなったモンテカルロは、モナコ公国の地区のひとつです。

（Google Map モンテカルロ）

モナコ自体はカジノが有名であり、確率を利用するところから
この名前がついたようです。
（wikipediaによるとジョン・フォン・ノイマンが命名したそうです。）

モンテカルロ法（Monte Carlo Method）は、
乱数を用いた数値計算手法の総称です。

例えば以前の記事で紹介した巡回セールスマン問題を、
「乱数でめぐる順番を生成→経路長が短くなったら更新」
を繰り返してとく場合、これもモンテカルロ法によるアプローチ
として扱うことができます。

2. 実装＆計算結果

円周率を求める仕組みのイメージをFig. 1に示します。

以下が実装したプログラムです。

import matplotlib.pyplot as plt
import random

図示したいのでmatplotlib、乱数使うのでrandomをimportします。

random.radom()が返してくれるのは0.0 <= x < 1.0のfloatです。
公式ドキュメント

ランダムな位置に打たれた点がx^2+y^2=1の内側かどうかを
7行目のif文で判断してカウントをしています。

N = 1000
pi_quarter = 0
plt.axes().set_aspect('equal')
for i in range(N):
    x = random.random()
    y = random.random()
    if (x**2 + y**2) < 1:
        pi_quarter += 1
        plt.scatter(x, y, c='blue')
    else:
        plt.scatter(x, y, c='green')
print(str(pi_quarter)+" / "+str(N))
print("predicted pi: "+str(pi_quarter*4/N))

plt.show()

上を実行すると↓

790 / 1000
predicted pi: 3.16

プログラム1行目のNの値を増やすとより正確な値が出やすくなります。

プログラムは以下からダウンロードできます。
Google Colab
Github

3. まとめ

• モンテカルロ法とは乱数を使った数値計算手法の総称
• 例題として円周率を求めてみた

以上です。

---