モンテカルロ法とは, 乱数計算, 円周率をpythonで計算

モンテカルロ法(Monte Carlo method)とは乱数を用いて数値計算を近似的に算出する方法のことです。例えばサイコロの1の目が出る確率は理論上は1/6ですが、それをシミュレーションで何万回もサイコロを振って、その出た目から確率を求めると1/6に非常に近くなる。というものです。本当にそうなるか具体例で説明します。

＜pythonによる実装具体例①＞
サイコロの目が出る確率として、1～6の値が出る確率を求めます。結果は以下のとおり、値が1/6に近い値になっております。

import numpy as np

a = np.random.randint(1,7,100000) # 1～6までの値を10万個出力する

print(np.count_nonzero(a==1)/100000) # 1の出る確率
print(np.count_nonzero(a==2)/100000) # 2の出る確率
print(np.count_nonzero(a==3)/100000) # 3の出る確率
print(np.count_nonzero(a==4)/100000) # 4の出る確率
print(np.count_nonzero(a==5)/100000) # 5の出る確率
print(np.count_nonzero(a==6)/100000) # 6の出る確率

⇒ 0.1653 　0.16735 　0.16653 　0.16817 　0.16664 　0.16601

＜pythonによる実装具体例②＞
モンテカルロ法を説明するのに有名な例です。以下の様に四角の中にランダムに点を打ち、円の中に入った点の数と、四角の中に入ってる点の数の比から円周率を求めるというものです。

実装プログラムは以下のとおり。for分を使わない方が高速に演算できますし、シンプルな記述になります。また円の4分の一で実施しても、同じ様に円周率を求めることができます。

import numpy as np
import matplotlib.pyplot as plt

sample = 100000                          # サンプル数
dot = np.random.rand(2,sample)    # ランダムサンプル
c_in = (dot[0]**2 + dot[1]**2) <= 1  # 円の中の点
c_out = np.logical_not(c_in)                # 円の外の点
count = np.count_nonzero(c_in)          # 円の中の点の数
pi = 4 * count/sample                        # 円周率
print(pi)

plt.scatter(dot[0][c_in], dot[1][c_in], s=0.1)
plt.scatter(dot[0][c_out], dot[1][c_out], s=0.1)
plt.show()

円周率は3.14176となり、円周率を近似的に求めることができました(結果はその都度変わります)。プロット結果は以下のとおり。

サブチャンネルあります。⇒ 何かのお役に立てればと