【Pythonで学ぶ】モンテカルロ法入門：平均・分散の計算

モンテカルロ法を、数値計算プログラミングの導入として学びたい
でもプログラミング苦手だし、どうやって勉強したらいいか分からない
手を動かしながらモンテカルロ法のプログラミングのコツを掴みたい

こういった悩みを抱えた読者さんを想定しています。
（私も学部生の頃はプログラミングが苦手であり、コツを掴むのに苦労していました笑）

この記事では、「モンテカルロ法」の導入的知識とプログラミングについて、Python言語を用いて説明します。

本記事を読み終える頃には
初学者の方でも、平均・分散の計算という初歩的な問題を解くことで、モンテカルロ法プログラミングの感覚を掴めるようになります。

本記事を執筆している、私トーマは、東京大学で物理学者として勤務し、理論物理や数値計算プログラミングに関する仕事をしています。

モンテカルロ法とは
実際に問題を解いてみよう
1. 考え方
2. 解答・ソースコード
まとめ
参考書籍

モンテカルロ法とは

モンテカルロ法とは「乱数を用いた統計的な計算手法」の総称です。
以下の節で、モンテカルロ法の乱数生成から、基本的な数値積分（モンテカルロ積分）について学んでいきます。

乱数生成

まず、モンテカルロ法に必要な乱数を生成していきます。
厳密にいうと、コンピュータは何らかの規則に沿った数列しか生成できないため、ある周期を持った疑似乱数が生成されることとなります。
（それでも、周期が十分長ければ、実用上は乱数と見做せます）

Pythonで乱数を生成するためには、以下のように、NumPyと呼ばれるライブラリパッケージを読み込む必要があります。

import numpy as np 	# numpy パッケージを np という名前で読み込む

そのパッケージ内のrandomモジュールを用いて、乱数を生成します。

一様乱数を生成したい場合

この場合は、以下のように書きます。

np.random.rand()       # [0,1)の範囲で一様乱数を生成する

標準正規分布に従う乱数を生成したい場合

この場合には、以下のように書きます。

np.random.randn()       # 標準正規分布に従う乱数を生成する

これらの乱数生成に使われる乱数ジェネレータはPCG-64というもので、\( 2^{128} \) という非常に長い周期をもちます。

モンテカルロ積分

以下のような積分
\begin{eqnarray}
I = \int_{a}^{b}dx \, f(x) \tag{1} \label{int}
\end{eqnarray}
をモンテカルロ法を用いて評価することを考えます。
（\( a < b \)である。）

ここで、ステップ関数 \( \Theta(x) \) を用いて
\begin{eqnarray}
P(x) \equiv \frac{\Theta(b-x) – \Theta(a-x)}{b-a} \tag{2} \label{prob}
\end{eqnarray}
と定義すると、\( P(x) \) は区間 \([a, b]\) での一様分布の確率密度関数と見做せます。

これを用いると、式 \eqref{int} は
\begin{eqnarray}
I = \int_{\infty}^{-\infty}dx \, P(x)(b-a) f(x) \equiv \int_{\infty}^{-\infty}dx \, P(x) \tilde{f}(x) \tag{3} \label{int2}
\end{eqnarray}
と書けます。
ただし、\( \tilde{f}(x) \equiv (b-a) f(x) \) としました。

式 \eqref{int2} を見てみると

\(I \) は \(\tilde{f}(x) \) の平均
分散 \(\sigma^{2}\) は \(\sigma^{2} = \displaystyle{\int^{\infty}_{-\infty}dx \, P(x) \left[ \tilde{f}(x) – I \right]^{2}} \)

となることが分かります。

一様乱数を用いた数値積分

区間 \([a, b]\) で一様に分布した（つまり \( P(x) \) に従う）\(N\) 個の乱数 \(r_{i}\) を用いて
\begin{eqnarray}
S = \frac{1}{N} \sum_{i=0}^{N-1} \tilde{f}(r_{i}) \tag{4} \label{sum}
\end{eqnarray}
と定義します。

この \(S\) は、母平均 \(I\) に対する標本平均と見做せるため、\(N\) が十分大きい極限で
\(S\) の分布は平均 \(I\) 分散 \(\sigma^{2}/N\) の正規分布に近づいていきます。
（これを中心極限定理という）

\(S\) を用いて \(I\) を近似的に求めようとした場合、かなりの高確率（約99.7%）で、誤差\( |I – S| \) は
\begin{eqnarray}
|I – S| < \frac{3\sigma}{\sqrt{N}} \tag{5} \label{Dif}
\end{eqnarray}
であることを意味しています。
\(S\) による積分評価の誤差は \(O \left( \displaystyle{\frac{\sigma}{\sqrt{N}}} \right)\) 程度であると言えます。

当然、\(N\) の値を大きくすれば、誤差は小さくなります。

実際に問題を解いてみよう

問題

一様乱数を用いたモンテカルロ積分により
\begin{eqnarray}
I = \int_{0}^{1}dx \, e^{-x} \tag{}
\end{eqnarray}
を評価せよ。さらに平均値と分散を計算し、理論値と比較せよ。

考え方

積分 \(I\) を母集団における平均とみなすと、\(P(x) = \displaystyle{\frac{\Theta(1-x) – \Theta(0-x)}{1-0}}\)、\(\tilde{f}(x) = (1-0)e^{-x}\)と書けるので

平均：\(I = 1 – e^{-1}\)
分散：\(\sigma^{2} = \displaystyle{\int_{-\infty}^{\infty} dx \, P(x) \left[\tilde{f}(x) – I \right]^{2} = -\frac{1}{2} \left(1 – e^{-1}\right) \left(1 – 3e^{-1}\right)} \)

つまり、\(N\) 個のサンプリング点を用いた評価値は、平均 \(I\) 分散 \(\sigma^{2}/N\) の正規分布に従うと見做せます。

たとえば、\(N=100\)サンプルの積分計算を \(M=1000\) 行い、その分布を図で確認してみましょう。

解答・ソースコード

分布の確認には、ヒストグラム（plt.hist 関数）を用います。
以下のように書けば、計算結果と解析解を重ねて表示できます。

import numpy as np 	# numpy パッケージを np という名前で読み込む
import matplotlib.pyplot as plt

N, M = 100, 1000        # N個のサンプリング点での積分をM回行う

X = np.random.rand(N, M)
ans = np.average(np.exp(-X), axis=0)    # M回分の積分結果を取得
ans = np.array(ans) - (1-np.exp(-1))    # 結果を母平均Iから測る
x = np.linspace(-0.1, 0.1, 200)
s2 = -0.5*(1-np.exp(-1))*(1-3*np.exp(-1))/N
y = np.exp(-x**2/s2/2)/np.sqrt(2*np.pi*s2)
plt.plot(x, y)     # 平均0、分散s2の正規分布をプロットする
plt.hist(ans, 51, (-0.1,0.1), density=True) 
             # density=Trueとすると、積分値が1となるように振幅が調整される
plt.show()

他にも、\(N, M\) の値を変えた時にどう評価値が変わるかもみてみると良いでしょう。