重点サンプリング (1)（2014/9/20）のつづき。今回は重点サンプリングを実装して、前回うまく計算できなかった積分が正しく計算できるか確認したい。

重点サンプリング

重点サンプリングの導出は、従来のモンテカルロ積分の式から簡単にできる。新しく重点関数 $g(x)$ を導入し、モンテカルロ積分の式を以下のように変形する。

$\displaystyle \mathbb{E}_{f} [ h(X) ] = \int_{\mathcal{X}} h(x) f(x) dx = \int_{\mathcal{X}} h(x) \frac{f(x)}{g(x)} g(x) dx = \mathbb{E}_{g} \biggl[ \frac{h(X) f(X)}{g(X)} \biggr]$

$g(x)$ が突然現れているが、分母と分子でキャンセルされるので追加しても問題ない。3式目から4式目の変形は、

$h(x) \frac{f(x)}{g(x)}$ が従来のモンテカルロ積分の式（上式の2式目）の $h(x)$ に対応する
$g(x)$ が従来のモンテカルロ積分の式（上式の2式目）の $f(x)$ に対応する

と考えると理解できる。重点サンプリングは、この結果を利用し、確率分布 g(x) から m 個のサンプル $(x_1, \cdots x_m)$ を生成し、上記の期待値の積分を次式で近似する。

$\displaystyle \bar{h}_{m} = \frac{1}{m} \sum_{j=1}^{m} \frac{f(x_j)}{g(x_j)} h(x_j)$

$m \to \infty$ のとき $\bar{h}_{m} \to \mathbb{E}_{f} [ h(X) ]$ となり、オリジナルの期待値の積分が求められる。重点サンプリングでは、オリジナルの $f(x)$ からサンプリングする代わりに重点関数 $g(x)$ でサンプリングしている。その代償として、 $h(x)$ を、重要度重み（importance weight） $r(x) = f(x) / g(x)$ で重み付けして補正しているという解釈ができる。重点サンプリングを実行するには、重点関数 $g(x)$ は $f(x)$ の代わりにサンプリングができる密度関数である必要がある。

重点サンプリングの例

というわけでさっそく重点サンプリングを実装してみよう。重点サンプリング (1)（2014/9/20）で通常のモンテカルロ積分ではうまく計算できなかった下の積分を重点サンプリングで計算してみる。

$\displaystyle \int_{-\infty}^{\infty} I(x \geq 5) N(x;0,1) dx$

重点関数 $g(x)$ として何を使うかが重要になりそうだ。とりあえず平均が5、分散が1の正規分布 $N(x;5,1)$ でやってみよう。

テキスト（例3.5）だと切断指数分布（truncated exponential distribution）を使えとか難しいこと書いてあるのだが、それは後でやってみよう。

下の図で青い正規分布がオリジナルの $f(x)$ で緑色の正規分布が重点関数 $g(x)$ である。

f:id:aidiary:20140921184557j:plain

実行結果は、

scipy.integrate: 2.86651570352e-07 2.86652756236e-07
normal monte carlo integration: 0.0
importance sampling: 2.96394456018e-07

普通のモンテカルロ積分では、積分が0になってしまうが、重点サンプリングを使うと2.96394e-07となり、scipy.integrateの結果と比較的近いことが確認できる。

青色の正規分布 $f(x)$ からサンプリングするとほとんどのサンプルが捨てられるけど、緑色の正規分布からサンプリングすると $h(x)$ とかぶっている半分のサンプルは無駄にならずに使われたんだろうなと推定できる。ただ、そのサンプルはあくまで $g(x)$ からのサンプルであるため、重要度重みで補正してあげる必要があるというわけか。

何となく重点関数は $h(x) f(x)$ の値が大きいところ（今回は5.0から6.0あたり）に高い密度がくるような分布を選べばよいのかなという推測がつく。次回は、重点関数をいろいろ変えたときにモンテカルロ積分の収束がどのように変わるか確認してみたい。