数学の具体的な計算にPythonを使って、数学もPythonも同時に学んでしまいましょう。今回はPythonを使った確率・統計の問題として、モンテカルロ法で円周率を求める問題をみてみたいと思います。本記事のPythonの学習事項としては、離散データの作成・処理の方法、散布図の作成を学びます。ただし、動作の確認はJupyterNotebookを用いた場合で行っています。インターフェイスなどの動作環境の違いによって適宜変更点があるかもしれません。
乱数
Pythonには（擬似）乱数が用意されています：
import random as rd rd.random() # 0.0以上1.0未満の小数をランダムに出力
そこで２つの乱数 を生成します。1.0 未満の２つの小数をランダムに出力したとき、そのデータ点を多数（その数を
としましょう）プロットしていくと、
という領域をまんべんなく埋め尽くしていくはずです。
実際、1,000個のデータ点をランダムにプロットするために
import random as rd def x(k): return rd.random() def y(k): return rd.random() P = 1000 R = [] for k in range(P): R.append([x(k),y(k)]) import matplotlib.pyplot as plt x, y = zip(*R) plt.scatter(x, y) plt.show()
と入力すると
のようになって の範囲にほぼ一様にデータ点が分布します。
モンテカルロ法
そこで、これらのデータ点のうち、 という条件を満たすデータ点の数（この数を
としましょう）を求めてみます。これは第１象限内の半径１の円のなかにあるデータ点の数に等しいので、データ点の数の比
をとれば、半径１の円の面積（つまり円周率です）の４分の１におおよそ等しいはずです：
P = 1000 import random as rd def x(k): return rd.random() def y(k): return rd.random() def c(k): if x(k)**2+y(k)**2 <= 1: return 1 else: return 0 s = [] for k in range(P): s.append(c(k)) S = sum(s) S/P*4 # 円周率と見比べるために４倍しておきました
もちろんデータ点が少なければ誤差が大きいわけですが、データ点を多くとり、試行回数を増やしていくほど良い近似値が得られるだろうと期待されます。これをPythonを使ってみていくことにしましょう。そこで以下では、10,000個のデータ点をとり、試行回数を500回に設定して、上記の範囲内に入ったデータ点の個数を考えてみます。
円周率の計算
P = 10000 Q = 500 import random as rd def x(k): return rd.random() def y(k): return rd.random() def c(k): if x(k)**2+y(k)**2 <= 1: return 1 else: return 0 T = [] for n in range(Q): s = [] for k in range(P): s.append(c(k)) def S(n): return sum(s) T.append(S(n))
T[n] は n 回目の試行で、10,000個のデータ点のうち、データ点が半径１の円の範囲内に入った場合の数です。
範囲内に入ったデータ点の個数ごとに、どれくらいの試行数でそれが実現したかをみるために、ヒストグラムを作成します：
import matplotlib.pyplot as plt plt.hist(T,bins=500, range=(7700, 8000)) plt.show()
すると
のようなヒストグラムが得られて、おおよそ7850のあたりにピークがあるように見えます。
平均値・標準偏差を求めると
import statistics as st m = st.mean(T) s = st.stdev(T) print('平均： {}'.format(m)) print('偏差： {}'.format(s))
より（いま乱数を使っているので、次の数値は再現されるものではありませんが）それぞれ
などと求まります。よって円周率の値が
と計算されます。結構良い精度で円周率が求められることがわかります。
今回は、乱数を使った数値計算の方法であるモンテカルロ法を使って円周率を計算しました。
