モンテカルロ積分の直感的理解と必要サンプル数の導出・精度向上などまとめ

乱数を用いたモンテカルロ法はモンテカルロ積分やMCMCなど、様々な応用に用いられます。モンテカルロ積分の計算手順自体はそれほど難しくない一方で、数式上は難しそうな表記がなされることが多いので当記事では直感的な理解ができるように取り扱いを行いました。
作成にあたっては「統計学実践ワークブック」の数式表記やWikipediaの定義などを参考に取りまとめを行いました。

日本統計学会公式認定統計検定準1級対応統計学実践ワークブック

3,080円(04/23 02:53時点)

Amazon

Contents

1 モンテカルロ積分の概要
2 モンテカルロ積分の精度
- 2.1 推定値の分散・標準偏差の計算
- 2.2 精度の向上
3 参考
- 3.1 関連

モンテカルロ積分の概要

モンテカルロ法とモンテカルロ積分

「統計学実践ワークブック」ではモンテカルロ法とモンテカルロ積分の使い分けがわかりにくいので、Wikipediaの内容の参照を行う。

In mathematics, Monte Carlo integration is a technique for numerical integration using random numbers. It is a particular Monte Carlo method that numerically computes a definite integral.
https://en.wikipedia.org/wiki/Monte_Carlo_integration

モンテカルロ積分(Monte Carlo integration)はモンテカルロ法(a particular Monte Carlo method)の一種とあり、「モンテカルロ法に内包される手法である」と考えると良い。

「モンテカルロ法」は「乱数を用いた計算」を指すので、「モンテカルロ積分」は「乱数を用いた積分値の計算」に対応すると理解しておくと良いと思われる。具体的な積分の式表示に関しては次項で取りまとめを行う。

$0 \leq x \leq 1$区間のモンテカルロ積分

$$
\large
\begin{align}
\theta = \int_{0}^{1} g(x) dx
\end{align}
$$

モンテカルロ積分は一様乱数を用いて上記の$\theta$の推定を行う手法である。ここでの積分の区間が$0 \leq x \leq 1$であることから、一様分布$U[0,1]$に基づいて乱数を生成し、これを用いて$\theta$を推定することを考える。

確率変数$X$が$X \sim U[0,1]$であるとき、$g(X)$の期待値$E[g(X)]$に関して下記が成立する。
$$
\large
\begin{align}
E[g(X)] &= \int_{0}^{1} g(x) f(x) dx \\
&= \int_{0}^{1} g(x) \times 1 dx = \int_{0}^{1} g(x) dx
\end{align}
$$

上記で用いた$f(x)$は一様分布$U[0,1]$の確率密度関数であるので、$0 \leq x \leq 1$の区間では$f(x)=1$であることを元に代入を行なった。ここで下記のように$g(X_1),…,g(X_m)$の標本平均$\hat{\theta}$を考える。
$$
\large
\begin{align}
\hat{\theta} = \frac{1}{m} \sum_{i=1}^{m} g(X_i)
\end{align}
$$

ここで$m \to \infty$のとき、大数の弱法則より$\displaystyle \frac{1}{m} \sum_{i=1}^{m} g(X_i) \to E[g(X)]$が成立する。これは$\displaystyle \lim_{m \to \infty} \hat{\theta} = \theta$と同義である。

ここまでは原理の確認を行なったが、実際にモンテカルロ積分を行う際は$U[0,1]$にしたがって$x_1,…,x_n$のサンプリングを行なったのちに、$g(x_1),…,g(x_n)$の標本平均を計算することで$\theta$の推定を行うと考えておくとよい。

下記にPythonを用いた$\displaystyle \theta = \int_{0}^{1} e^x dx$の推定の例を示す。

import numpy as np

np.random.seed(20)

x_100 = np.random.rand(100)
x_1000 = np.random.rand(1000)
x_10000 = np.random.rand(10000)
x_100000 = np.random.rand(100000)
y_100 = np.e**x_100
y_1000 = np.e**x_1000
y_10000 = np.e**x_10000
y_100000 = np.e**x_100000


print("Population theta: {:.5f}".format(np.e-1))
print("n=100, Estimated theta: {:.5f}".format(np.mean(y_100)))
print("n=1000, Estimated theta: {:.5f}".format(np.mean(y_1000)))
print("n=10000, Estimated theta: {:.5f}".format(np.mean(y_10000)))
print("n=100000, Estimated theta: {:.5f}".format(np.mean(y_100000)))

・実行結果

Population theta: 1.71828
n=100, Estimated theta: 1.73545
n=1000, Estimated theta: 1.70142
n=10000, Estimated theta: 1.73034
n=100000, Estimated theta: 1.71785

実行結果より$n=100,000$の精度が一番高いことが確認できる。また、$n=10000$の精度はあまり良くないが、乱数の設定によっては十分生じうると考えておくと良いと思われる。

区間が$0 \leq x \leq 1$でない場合のモンテカルロ積分

前項では$0 \leq x \leq 1$区間におけるモンテカルロ積分に関して確認を行なったが、当項では区間を変えた際のモンテカルロ積分をどのように計算するかについて確認を行う。

$$
\large
\begin{align}
\theta = \int_{0}^{3} e^{x} dx
\end{align}
$$
以下、具体的に上記のように$g(x)=e^{x}$の$0 \leq x \leq 3$の区間での積分値の推定に関して考える。

$$
\large
\begin{align}
\int_{0}^{3} e^{x} dx = \int_{0}^{1} e^{x} dx + \int_{1}^{2} e^{x} dx + \int_{2}^{3} e^{x} dx
\end{align}
$$
積分値の推定にあたっては上記のように積分の式を分解して考えることで、モンテカルロ積分を行う際に用いる確率密度関数を該当区間で$f(x)=1$となる一様分布に設定することができ、そのまま推定値を計算するだけで良いのでシンプルに考えることができる。

前項のPythonの処理に基づいて$\displaystyle \theta = \int_{0}^{3} e^x dx$の推定の例を示す。

import numpy as np

np.random.seed(20)

theta = 0
for i in range(3):
    x = np.random.rand(100000) + i
    y = np.e**x
    theta += np.mean(y)


print("Population theta: {:.5f}".format(np.e**3-1))
print("Estimated theta: {:.5f}".format(theta))

・実行結果

Population theta: 19.08554
Estimated theta: 19.09739

モンテカルロ積分の精度

推定値の分散・標準偏差の計算

$$
\large
\begin{align}
\hat{\theta} = \frac{1}{m} \sum_{i=1}^{m} g(X_i)
\end{align}
$$

モンテカルロ積分は上記の式を元にパラメータ$\hat{\theta}$の推定値を計算する手法である。これは$\theta$は$\displaystyle \lim_{m \to \infty} \hat{\theta} = \theta = E[g(X_i)] = E[g(X)]$のように考えることができることに基づくが、ここで推定量・推定値に対応する$\hat{\theta}$に対して、$\hat{\theta}$の分散$V[\hat{\theta}]$や標準偏差$SD[\hat{\theta}]$は推定値の精度に対応すると考えることができることは抑えておくと良い。

以下$g(x) = e^x$に対して$X_1,…,X_m \sim U[0,1]$を用いて$\displaystyle \theta = \int_{0}^{1} e^x dx$を推定する際の標準偏差に関して考える。先に$g(X)$の分散$V[g(X)]$の計算を行う。
$$
\large
\begin{align}
V[g(X)] &= E[(g(X)-E[g(x)])^2] \\
&= \int_{0}^{1} \left( g(x) – \int_{0}^{1} g(x) dx \right)^2 dx \\
&= \int_{0}^{1} \left( e^x – \int_{0}^{1} e^x dx \right)^2 dx \\
&= \int_{0}^{1} \left( e^x – \left[ e^x \right]_{0}^{1} \right)^2 dx \\
&= \int_{0}^{1} ( e^x – e + 1 )^2 dx \\
&= \sigma_g^2
\end{align}
$$
上記では結果を$\sigma_g^2$のようにおいたが、推定値の標準偏差を考える場合は実際には値が与えられる場合や推定を行なった値を用いると考えておくとよい。

次に$X_1,…,X_n$が$i.i.d.$のとき、分散$V[\hat{\theta}]$を上記の$\sigma_g^2$を用いて表すことを考える。
$$
\large
\begin{align}
V[\hat{\theta}] &= V \left[ \frac{1}{m} \sum_{i=1}^{m} g(X_i) \right] \\
&= \frac{1}{m^2} V \left[ \sum_{i=1}^{m} g(X_i) \right] \\
&= \frac{1}{m^2} \times m \times V[g(X_i)] \\
&= \frac{1}{m} V[g(X_i)] \\
&= \frac{\sigma_g^2}{m}
\end{align}
$$

上記より、推定量・推定値の標準偏差$SD[\hat{\theta}]$は下記のように表すことができる。
$$
\large
\begin{align}
SD[\hat{\theta}] &= \sqrt{V[\hat{\theta}]} \\
&= \frac{\sigma_g}{\sqrt{m}}
\end{align}
$$

・参考
期待値、分散の公式

精度の向上

参考

・自然科学の統計学解答
 ・統計学実践ワークブック解答