モンテカルロ法による積分（１） - 機械学習に詳しくなりたいブログ

これからマルコフ連鎖モンテカルロ法を勉強していこうと思っています。今回はその第一歩。モンテカルロ法というのは何か特定の手法を表すのではなく、乱数を用いて数値計算すること全般を指すようです。ですのでモンテカルロ法による積分とは、乱数を用いて積分を数値計算する、ということです。

$f(x)$ が確率分布 $p(x)$ に従うなら、期待値は

$\displaystyle E[f(x)] = \int f(x)p(x)dx \tag{1}$

と書けます。いま、ある関数 $f(x)$ の積分 $\displaystyle \int_{a}^{b} f(x)dx$ を考えます。区間 $[a,b$ ]で一様分布する確率密度関数 $p(x)$ は、

$\displaystyle p(x)= \begin{cases} \displaystyle \frac{1}{b-a} & (a \le x \le b) \\ 0 & (otherwise) \end{cases}\tag{2}$

ですから、これを使って

$\begin{eqnarray*} \displaystyle \int_{a}^{b}f(x)dx &=& \int_{a}^{b}\frac{f(x)}{p(x)}p(x)dx \tag{3} \\ &=&(b-a) \int_{a}^{b}f(x)p(x)dx \tag{4}\\ &=&(b-a) E[f(x)] \tag{5} \end{eqnarray*}$

と変形できます。 $p(x)$ は $x$ に依存しませんから積分の外に出せます。積分から一部を出して一部を残すっていう不思議な式変形でちょっと戸惑いそうです。

さて、大数の法則*1より、 $N$ が十分大きければサンプル平均は期待値に近づいていきますから、

$\displaystyle \int_{a}^{b} f(x)dx \simeq \frac{b-a}{N}\sum_{n=1}^{N}f(x^{(n)}) \tag{6}$

と近似できます。ここで $x^{(n)}$ は $p(x)$ からサンプルした値です。よって、積分区間で一様分布する $p(x)$ からサンプル $x^{(n)}$ をとることができれば、 $f(x)$ の積分が期待値計算により求められることになります。この考え方で積分の数値計算をする方法をモンテカルロ積分と呼びます。最初に知ったときはなるほどなあと随分感心してしまいました。

さて、モンテカルロ積分で

$\displaystyle \int_{0}^{4}\sqrt{2x+1}dx \tag{7}$

を求めてみます。以下の図のような関数、積分区間です。 f:id:opabinia2:20180303010504p:plain 区間 $[0,4$ ]で一様分布する確率密度関数 $\displaystyle p(x)=\frac{1}{4}$ を用いて、

$\begin{eqnarray*} \displaystyle \int_{0}^{4}\sqrt{2x+1}dx &=& \int_{0}^{4}\frac{\sqrt{2x+1}}{p(x)}p(x)dx\tag{8}\\ &=&4\int_{0}^{4}\sqrt{2x+1}p(x)dx\tag{9}\\ &\simeq& \frac{4}{N}\sum_{n=1}^{N} \sqrt{2x^{(n)}+1} \tag{10} \end{eqnarray*}$

です。 $N$ の値を変化させて、式(7)の解： $\displaystyle \frac{26}{3}$ との誤差を確認してみました。 f:id:opabinia2:20180303014847p:plain

きちんと近似できていることが確認できました。サンプル数 $N$ が2桁増えると誤差が1桁減るような傾向です。

今回のコードです。サンプル数が少ないと値が安定しなかったので、繰り返し計算して平均を求めています。対数グラフで綺麗な直線が見たかったので繰り返し数は多め。

# モンテカルロ積分テスト

import numpy as np
import matplotlib.pyplot as plt


# 積分対象の関数
def f(x):
    return np.sqrt(2*x+1)


# Nサンプル発生させて積分を近似計算し、真値との誤差を計算（積分の真値は26.0/3.0）
def err(N):
    p = np.random.uniform(0, 4, N)
    return np.abs(4*np.average(f(p)) - 26.0/3.0)


# ランダムシードを固定
np.random.seed(0)

# M回誤差を計算して平均をとる
M = 100

error = []

# サンプリング数のリスト
N = [10**i for i in range(2, 9)]

for i in N:
    # 計算時間調整
    # if i >= 1e5:
    #     M = 10
    error.append(np.average([err(i) for j in range(M)]))

plt.yscale("log")
plt.xscale("log")
plt.grid(which="both")
plt.xlabel("N")
plt.ylabel("error")
plt.plot(N, error)
plt.show()

*1:大数の法則は弱法則、強法則とあってややこしい。高校の数学で習う程度の理解で困ることはなさそうか？いつかきちんと理解したい