Pythonによるモンテカルロ法（2014/6/20）のつづき。今回は、重点サンプリング（importance sampling）の実験をしてみたい。重点サンプリングは、PRMLの11章にも出ている。重点サンプリングを使うと期待値のモンテカルロ積分をより効率的に（＝少ないサンプル数で高精度に）計算できる。

記事が長くなりそうなので今回はイントロのみということでなぜ重点サンプリングが必要かという前置きだけ書き記しておきたい。

普通のモンテカルロ積分

モンテカルロ積分（2014/7/28）から続く何本かの記事でまとめたけど、下記の期待値の積分をサンプリングによって求めるのがここでの目的。

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

ここで、 $f(x)$ は確率分布、 $h(x)$ は任意の関数である。モンテカルロ積分では、確率分布fからm個のサンプル $(x_1, \cdots x_m)$ を生成し、上記の期待値の積分を次式で近似する。

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

このアルゴリズムは、サンプルを確率分布 $f(x)$ から生成できるというのが前提になっている。しかし、 $f(x)$ からサンプリングできなかったり、 $f(x)$ からのサンプルを使って近似値を計算すると非常に効率が悪い（＝大量のサンプルを生成しないと近似精度が上がらない）ケースがあるというわけ。

効率が非常に悪いケース

たとえば、下の積分を通常のモンテカルロ積分で計算しようとすると非常に効率が悪い。

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

ここで、 $N(x;0,1)$ は標準正規分布。上の式を期待値の積分の式に合わせて書き直すと

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

となる。ここで、I([条件]) は、条件を満たす範囲で1を返し、それ以外で0を返す関数とする。この期待値の積分を通常のモンテカルロ法で求めてみよう。

実行結果は、

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

となる。正しい答えは、2.8665e-07なのに0.0となってしまった。まあほとんど合ってるんだけど（笑）正確には違う。それは図示してみるとよくわかる。

f:id:aidiary:20140920190003p:plain

図示してみると、一見、 $f(x)$ と $h(x)$ はほとんど重なりがないので被積分関数 $h(x)f(x)$ の全区間の積分は0でいいんでない？と思うかもしれないけど、図の丸で囲んだ部分を拡大してみるとちゃんと面積があるのだ！つまり、0ではない。

それでは、通常のモンテカルロ積分で積分が0になってしまう理由を考察してみよう。モンテカルロ積分では、 $f(x)$ （下の図の青色の正規分布）からサンプル $x$ を生成し、そのサンプル $x$ を $h(x)$ の関数に通してから下式で和を取って積分の近似値を求める。

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

f:id:aidiary:20140920190006j:plain

この際、上の図の黒い矢印の範囲から生成されたサンプル $x$ は、 $h(x)=0$ となってしまうため最終的に積分の近似値の和の計算では使われないことになる。結局、積分の近似値の和の計算で実際に意味があるのは、 $h(x)$ が0とならない上の図の赤い矢印の範囲のサンプルだけというわけ。

しかし、赤い矢印の範囲のサンプルが生成される確率はとてつもなく低い。なぜかというとサンプルは青色の正規分布に従って生成されるため正規分布の裾野の裾野のまた裾野くらいの赤い範囲のサンプルが生成される確率はほとんど0なのだ。よって、サンプルスクリプトのように $N=1000$ サンプル程度生成しただけでは、赤い矢印の範囲のサンプルが生成される可能性はほぼなく、最終的な積分の近似値も0になってしまうというわけ。まあ、Nをものすごーく大きくすれば運良く赤い範囲のサンプルがいくつか生成されて、積分を近似できるようになるかもしれないけど効率はものすごーく悪くなってしまう。

PRMLの重点サンプリング 11.1.4 の節には次のような記述があった。この問題は高次元空間だとより深刻になるようだ。

我々の興味があるような確率分布では、しばしば $\mathbf{z}$ 空間の比較的小さい領域に大きな確率が集中している。この場合、高次元の問題では、非常に小さな割合のサンプルだけが和に対して大きな寄与をするため、一様サンプリングはとても非効率的である。 $p(\mathbf{z})$ が大きな領域、理想的には $p(\mathbf{z}) f(\mathbf{z})$ が大きな領域に入るサンプル点を本当は抽出したいのである。

（注）Rによるモンテカルロ法入門と記号の使い方が違うので混乱しやすい。PRMLでは、確率分布がp（本ブログではf）で任意の関数がf（本ブログではh）になっている。