モンテカルロ積分 (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

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

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