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

モンテカルロ法による積分（２）無限区間の積分で、正規分布の確率密度関数の積分をしました。正規分布からサンプリングしてモンテカルロ積分をしましたが、対象の積分区間 $[2,\infty$ ]に対して、正規分布全体からサンプリングするのは実はとても効率が悪い方法のようです。正規分布は $\pm 2\sigma$ の範囲に約95%が集中していますから、標準正規分布で2以上がサンプリングできる確率は約2.5%程度しかありません。つまりほとんどのサンプルは計算に必要のない無駄なものとなっていました。この問題の解決するのが重点サンプリングです。以下の式変形のように、元の確率分布の効率が悪ければ、効率の良い任意の確率分布 $q(x)$ を導入し、そこからサンプリングすれば良いわけです。この $q(x)$ は提案分布と呼ばれます。

$\begin{eqnarray*} \displaystyle E[f(x)] &=& \int f(x)p(x)dx \tag{1}\\ &=& \int \frac{f(x)p(x)}{q(x)}q(x)dx \tag{2}\\ &\simeq& \frac{1}{N}\sum_{n=1}^{N} \frac{f(x^{(n)})p(x^{(n)})}{q(x^{(n)})} \tag{3} \end{eqnarray*}$

この式変形はシンプルですけど、なるほどなあという感じですよね。

さて、先回は正規分布 $f(x)$ からサンプリングし、

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

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

として計算しましたが、提案分布を用いて重点サンプリングしてみます。提案分布は、以下のグラフのような2種類を使ってみました。（赤のプロットは積分対象で、青と緑のプロットが提案分布です）

f:id:opabinia2:20180305213650p:plain

1つは $\mu=2$ の正規分布。中心値が積分範囲の下限なので、これで50％のサンプルが有効に使われることになります。もう1つは指数分布を $x$ 軸に2だけシフトしたもの。これなら100％のサンプルが使われます。それぞれサンプル数 $N$ を変化させて真値との誤差を計算してみた結果が次のグラフです。

f:id:opabinia2:20180305213704p:plain

指数分布を用いた場合、同じ精度を得るのに約1/100のサンプル数で済みました。 $\mu=2$ の正規分布がその中間くらい。 $\mu=3$ とかにして、50％よりもっと多くのサンプルが使われるようにしたほうが良いかと思って試してみましたが、 $\mu=2$ より僅かに効率が悪くなりました。サンプルの効率も重要ですが、積分対象の形に近いこと、というのも重要な要素のようです。また別の例で実験してみたいと思います。でもグラフで図示できれば確認できますが、高次元の場合は提案分布がその計算に適しているのか？判断するのは難しそうですね。

今回でモンテカルロ積分について書くのは3回目ですが、読み返してみると、その時々で例えば $p(x)$ が確率分布だったり、積分対象だったり、文字の使い方がわかりにくいですね、、、。

# 重点サンプリング 正規分布の確率密度関数の積分
import numpy as np
import matplotlib.pyplot as plt


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


# 標準正規分布確率密度関数
def normal(x):
    return (1/np.sqrt(2*np.pi))*np.exp(-(x**2)/2)


# λ=1 、2シフトの指数分布確率密度関数
def exponential(x):
    return np.exp(-(x-2))


# Nサンプル発生させて誤差を計算（積分の真値は約0.0227501）
def err(N):
    x = np.random.exponential(1, N) + 2
    # x = np.random.normal(2, 1, N)
    ret = vec_f(x)*(normal(x)/exponential(x)) # normal(x-2))
    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, 7)]

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.show()