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

モンテカルロ法による積分（１）の続きです。

今回は標準正規分布 $f(x)=N(0,1)$ の積分を考えます。 $\displaystyle \int_{a}^{b} f(x)dx$ であれば、先回と同じく一様分布 $p(x)$ を用いて $\displaystyle (b-a) \int_{a}^{b}f(x)p(x)dx$ として $p(x)$ からサンプリングすれば良いです。ただし、積分区間に無限を含む $\displaystyle \int_{2}^{\infty} f(x)dx$ の場合は同じように計算できません。無限区間で一様分布する確率密度関数は存在しませんので。

f:id:opabinia2:20180304004800p:plain

今回の場合は積分対象が正規分布ですから、 $f(x)$ から直接サンプリングすることが可能です。 $2 \le x$ のときにのみ1となる関数 $p(x)$ を使えば、

$\begin{eqnarray*} \displaystyle \int_{2}^{\infty} f(x)dx &=& \int_{-\infty}^{\infty} f(x)p(x)dx \tag{1}\\ &\simeq& \frac{1}{N}\sum_{n=1}^{N} p(x^{(n)}) \tag{2} \end{eqnarray*}$

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

です。よって $f(x)$ からサンプリングし、 $p(x)$ の平均を求めれば良いことになります。サンプル数Nを変化させて真値との誤差を計算してみました。 f:id:opabinia2:20180304002641p:plain

ちゃんと計算できていますね。適当な区間で打ち切って一様分布からサンプリングしてもそれなりの精度になりそうですけど。誤差の桁は先回のモンテカルロ法による積分（１）の実験より小さいですが、真値の桁数に依存するものです。

なお積分対象からサンプリングできなければ、式(3)の $p(x)$ と、適当な分布 $g(x)$ を用意して

$\begin{eqnarray*} \displaystyle \int_{2}^{\infty} f(x)dx &=& \int_{-\infty}^{\infty} \frac{f(x)p(x)}{g(x)}g(x)dx \tag{4}\\ &\simeq& \frac{1}{N}\sum_{n=1}^{N} \frac{f(x^{(n)})p(x^{(n)})}{g(x^{(n)})} \tag{5} \end{eqnarray*}$

とすればいいだけです。（参考：重点サンプリング）　積分区間でサンプリングできればどんな分布でも一応計算できるはずですが、積分対象とあまりにもかけ離れた形状の分布だと精度が出ないと思います。後日検証。

今回のコードです。今回も対数グラフで綺麗な直線が見たかったのでサンプル数は過剰です。一様分布より正規分布のほうが処理が遅いんでしょうか？先回より計算時間が相当かかりました。寝て起きたら終わってたので正確にはわかりませんが、多分5，6時間とか。サンプル数10E8を100回繰り返すとか、やりすぎでした、、、。

# モンテカルロ積分テスト 正規分布の確率密度関数の積分
import numpy as np
import matplotlib.pyplot as plt


# 積分対象の関数
def f(x):
    if x > 2:
        return 1
    else:
        return 0


# Nサンプル発生させて誤差を計算（積分の真値は約0.0227501）
def err(N):
    ret = vec_f(np.random.normal(0, 1, N))
    return np.abs(np.average(ret) - 0.0227501)


vec_f = np.vectorize(f)

# ランダムシードを固定
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()