pianofisica

Mathematics & Physics, Maxima, a bit Python & Wolfram, and Arts

Pythonで学ぶ確率・統計(モンテカルロ法)

数学の具体的な計算にPythonを使って、数学もPythonも同時に学んでしまいましょう。今回はPythonを使った確率・統計の問題として、モンテカルロ法で円周率を求める問題をみてみたいと思います。本記事のPythonの学習事項としては、離散データの作成・処理の方法、散布図の作成を学びます。ただし、動作の確認はJupyterNotebookを用いた場合で行っています。インターフェイスなどの動作環境の違いによって適宜変更点があるかもしれません。



乱数

Pythonには(擬似)乱数が用意されています:

import random as rd
rd.random()    # 0.0以上1.0未満の小数をランダムに出力

そこで2つの乱数  (x,y) を生成します。1.0 未満の2つの小数をランダムに出力したとき、そのデータ点を多数(その数を  P としましょう)プロットしていくと、 0\leq x\leq1, 0\leq y\leq1 という領域をまんべんなく埋め尽くしていくはずです。

実際、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()

と入力すると

f:id:pianofisica:20200904175544p:plain

のようになって  0\leq x\leq1, 0\leq y\leq1 の範囲にほぼ一様にデータ点が分布します。



モンテカルロ法

そこで、これらのデータ点のうち、  x^2+y^2\leq1 という条件を満たすデータ点の数(この数を  S としましょう)を求めてみます。これは第1象限内の半径1の円のなかにあるデータ点の数に等しいので、データ点の数の比  S/P をとれば、半径1の円の面積(つまり円周率です)の4分の1におおよそ等しいはずです:

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  # 円周率と見比べるために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個のデータ点のうち、データ点が半径1の円の範囲内に入った場合の数です。

範囲内に入ったデータ点の個数ごとに、どれくらいの試行数でそれが実現したかをみるために、ヒストグラムを作成します:

import matplotlib.pyplot as plt
plt.hist(T,bins=500, range=(7700, 8000))
plt.show()

すると

f:id:pianofisica:20200916133927p:plain

のようなヒストグラムが得られて、おおよそ7850のあたりにピークがあるように見えます。




平均値・標準偏差を求めると

import statistics as st
m = st.mean(T)
s = st.stdev(T)
print('平均: {}'.format(m))
print('偏差: {}'.format(s))

より(いま乱数を使っているので、次の数値は再現されるものではありませんが)それぞれ

\quad\displaystyle{\mu=7853.5 \qquad\quad \sigma=39.9}

などと求まります。よって円周率の値が

\quad\displaystyle{\pi=\frac{7853.5\pm39.9}{10000}\times 4=3.1414\pm0.0160}

と計算されます。結構良い精度で円周率が求められることがわかります。


今回は、乱数を使った数値計算の方法であるモンテカルロ法を使って円周率を計算しました。



キーワードMaxima、確率・統計、乱数、モンテカルロ法、データ作成、プロット