2014-09-21

重点サンプリング (2)

機械学習 PRML

重点サンプリング (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あたり）に高い密度がくるような分布を選べばよいのかなという推測がつく。次回は、重点関数をいろいろ変えたときにモンテカルロ積分の収束がどのように変わるか確認してみたい。

参考

2014-09-20

重点サンプリング (1)

機械学習 PRML

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）になっている。

先の例に当てはめると被積分関数 $h(x)f(x)$ が他と比べて大きな値をとる5.0から6.0くらいの範囲のサンプルを本当に抽出したいのだが、 $f(x)$ が標準正規分布だとそれができない！というのが問題になる。

重点サンプリングは、この問題を解決するため欲しい範囲のサンプルを重点的にサンプリングできる重点関数というものを導入しましょうというアイデア。重点サンプリングを使うと上の積分も0にならず、真の値に近い与えが得られる。

長くなったので重点サンプリングの具体的な実装はまた次回！

2014-08-30

モンテカルロ積分の収束テスト

機械学習 PRML

モンテカルロ積分 (2)（2014/8/20）のつづき。

前回の正規・コーシー-ベイズ推定量の積分のサンプリング方法には2通りの方法があった。コーシー分布からサンプリングする場合と正規分布からサンプリングする場合だ。どちらでもサンプリング数を多くすればほぼ同じ値に収束していたが、どちらがより効率的かを詳しく検証をしてみた。

前回までは、PRMLの表記に合わせて式を書いていたが、『Rによるモンテカルロ法入門』との表記の違いで混乱しやすいため今回からは後者に合わせていく。

下記の期待値の積分をサンプリングによって求めるのが目的。

$\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)$

$m \to \infty$ のとき $\bar{h}_{m} \to \mathbb{E}_{f} [ h(X) ]$ となる。収束判定ではよく使う手であるが、横軸にサンプル数を取り、縦軸に推定値 $\bar{h}_{m}$ をプロットして本当に $\mathbb{E}_{f} [ h(X) ]$ に収束するか確かめればよいだろう。あとサンプルから求まる以下の分散（スクリプトではルートをとった標準偏差）もついでにプロットしてみよう。

$\displaystyle v_m = \frac{1}{m^{2}} \sum_{j=1}^{m} [ h(x_j) - \bar{h}_m ]^{2}$

例3.3の積分の収束テスト

f:id:aidiary:20140830202932p:plain

実装の工夫だが、全サンプルを生成した後に各 m サンプルまでの和をcumsum()で求めている。このようなVectorizationの演算はforループで書くよりずっと高速だが慣れないと何をやっているのかわかりづらい。結果のグラフを見るとサンプル数を増やすことで期待値の積分の値に収束していくことが確認できる。

練習問題3.1の積分の収束テスト

f:id:aidiary:20140830202934p:plain

こっちは、冒頭で述べた正規・コーシー-ベイズ推定量の積分の場合。上のグラフがコーシー分布からサンプリングした場合で下のグラフが正規分布からサンプリングした場合。この結果からわかるように正規分布からサンプリングした方が収束が早いし、分散も小さいことがわかる。コーシー分布は正規分布に比べて裾野が広いため平均から遠く離れた極端なサンプルが生成されやすいことが収束の遅さにつながっているのだと思う。

モンテカルロ法ではどの分布からサンプリングするかも重要な要素なわけね。これが次の重点サンプリングにつながっていくのかな？

参考

Rによるモンテカルロ法入門読書ノート：3章モンテカルロ積分その2- Rによる実装

人工知能に関する断創録

このブログでは人工知能のさまざまな分野について調査したことをまとめています（更新停止: 2019年12月31日）

重点サンプリング (2)

重点サンプリング

重点サンプリングの例

参考

重点サンプリング (1)

普通のモンテカルロ積分

効率が非常に悪いケース

モンテカルロ積分の収束テスト

例3.3の積分の収束テスト

練習問題3.1の積分の収束テスト

参考