重点サンプリング (4)

重点サンプリング (3)（2014/10/2）のつづき。今回は、Rによるモンテカルロ法入門のテキストの例題をPythonで解いてみたい。まずは、練習問題3.4から。

fを標準正規分布 とし、 を として期待値 を計算せよ。

という問題。

• 標準正規分布 からサンプリングする通常のモンテカルロ積分
• 重点関数として一様分布 を使ったモンテカルロ積分

の収束性を比較してみたい。

(注) テキストでは、重点関数として一様分布 を使えと書いてあるが、この設定だとあとで見るようにサンプリングの効率が極端に悪く、正しい積分が得られないようだ。

青色が標準正規分布 、緑色が任意分布の 、赤色が被積分関数 、水色が重点関数 となっている。被積分関数のところを重点的に覆えるように重点関数を選んでみた。今まで (1)～(3) の記事を通して重点サンプリングの手順はさんざん書いてきたので今回はスクリプトと結果だけ。

実行結果は、

scipy.integrate: 0.0746157703288
normal monte carlo integrationi: 0.0731957831939
importance sampling: 0.0732103112973

どちらも同じくらいの値に収束している。ただ、ランダム性が入ってくるため実行するたびに結果は変わってしまう。収束性のグラフも表示してみよう。

この結果を見てわかるように重点サンプリングを使った方が収束が早いことが確認できる。ただ、この例題の結果を見ていて気になったけど、サンプル数を10万くらいまで多くしても結局のところモンテカルロ積分でscipy.integrateの値に有効数字3桁で一致する結果がなかなか得られなかった。こんなものなのかな？

ちなみにテキストに書かれているように一様分布 で実行した結果は下のようになった。重点サンプリングでは、一様分布から得られたサンプルxを用いて計算したの結果がほぼ0になるため積分の値も0になってしまう。何か勘違いしているかな？テキストが誤りなのかな？ご存知の方がいたら教えてください。