重点サンプリング (5)
重点サンプリング (4)(2014/10/11)のつづき。今回は、Rによるモンテカルロ法入門のテキストの例3.5を解いてみたい。
の積分を重点サンプリングで計算しろというもの。似たような積分は (1)~(3) までですでに取り上げてきた。そこでは、標準正規分布からのサンプリングでは4.5以上のサンプルを得ることはほぼありえないため、重点関数として平均値をずらした正規分布を用いてサンプリングした。今回は、重点関数として平均値をずらした正規分布ではなく、4.5で切り詰めた指数分布を使えというのが課題になる。
切断指数分布
切断指数分布(truncated exponential distribution)という分布は初めて聴いたのだが、指数分布から4.5以上の区間を取り出してきて、その区間の面積が1になるように再調整した分布ということがわかった。厳密な定義は、WikipediaのTruncated distributionの項を参照のこと。指数分布を4.5で切り詰めた分布は、
で表される。Wikipediaの定義からもちゃんと導出できることは確認した。条件として4.5より小さい範囲では、 となるので注意。早速、どのような分布なのか確認してみよう。
def g(x): if x < 4.5: return 0 return np.exp(-(x - 4.5)) ix = np.arange(-5, 15, 0.01) iy = [g(x) for x in ix] plt.plot(ix, iy, label="g(x)") plt.legend(loc="best") plt.show() print scipy.integrate.quad(g, -np.inf, np.inf)[0]
全区間で積分すると1になることも確認できた。
切断指数分布を用いた重点サンプリング
次に重点関数として4.5で切り詰めた切断指数分布を用いてモンテカルロ積分を計算してみよう。ついでに重点サンプリング (2)(2014/9/21)で用いた正規分布の結果もいっしょに収束性の違いを評価してみた。青色が本来のサンプリング分布である標準正規分布。緑色が1つめの重点関数候補である平均をずらした正規分布。赤色が2つめの重点関数候補である切断指数分布。
切断指数分布の関数は、引数に numpy.array を取れるようにするため苦しい実装をしている。forループで回すのは遅いためもっとよい実装があるかもしれない。また、切断指数分布は自分で実装した関数であるためscipyの標準分布のように乱数を生成できない。そのため、標準指数分布のサンプルに4.5を足して右にずらすことで切断指数分布のサンプルを生成している。どちらも同じ結果になる。
(注)scipyにも切断指数分布の関数があった scipy.stats.truncexpon。関数の定義を見ると切断は[0,b]であって、始点側での切断位置を指定できないみたいだ。
scipy.integrate: 3.39714902293e-06 normal monte carlo integration: 0.0 importance sampling (norm): 3.4405273784e-06 importance sampling (truncated expon): 3.32877254612e-06
通常のモンテカルロ積分では、4.5以上のサンプルxが正規分布から生成できないためほとんどのサンプルで となってしまい積分値が0になってしまう。重点サンプリングを使うとその問題がなくなる。
重点サンプリングで、重点関数に正規分布を使った場合と切断指数分布を使った場合を比較すると、切断指数分布を使った方が収束性が早く、分散も小さいことがわかる。おそらく、正規分布では4.5より小さい左半分のサンプルが無駄になるのに対し、切断指数分布では必ず4.5以上のサンプルが生成され、無駄になるサンプルがないためだと考えられる。
重点サンプリング (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になってしまう。何か勘違いしているかな?テキストが誤りなのかな?ご存知の方がいたら教えてください。
重点サンプリング (3)
重点サンプリング (2)(2014/9/21)のつづき。前回からけっこう日があいてしまったが、予定どおり重点関数をいろいろ変えたときにモンテカルロ積分の収束性がどのように変わるか検証した。 対象とする積分は下の式。通常のモンテカルロ積分では、平均0、標準偏差1の標準正規分布 からサンプリングするがそれだとほとんどのサンプルが使われなくなり積分が0になってしまう点が問題となっていた。
重点サンプリングでは、オリジナルの標準正規分布 の代わりに重点関数 からサンプリングして以下の値で積分を近似する。
以下の6つの正規分布を重点関数 として収束速度を調べてみた。
- 平均0、標準偏差1の正規分布(オリジナルの分布と同じ)
- 平均5、標準偏差1の正規分布
- 平均10、標準偏差1の正規分布
- 平均0、標準偏差3の正規分布
- 平均5、標準偏差3の正規分布
- 平均10、標準偏差3の正規分布
結果は、下のようなグラフになる。青色は scipy.integrate で求めた値で理論値にかなり近い積分値である。
この結果を考察してみると・・・
平均0、標準偏差1の場合は、ほとんどのサンプルが5以上の値を取ることがなく となってしまい、いくらサンプリングしても積分が0のままであることがわかる。前回、実験したように被積分関数にほどよく重なる平均5、標準偏差1のようにすると非常に早く理想値に収束する。ただし、平均10、標準偏差1の結果を見るとわかるように被積分関数と重なるからと言って、オリジナルの分布 から遠く外れていると、今度は となってしまい積分が0のままになる。下の段のように標準偏差を少し広げて、被積分関数とオリジナルの関数の両方に少しずつ重なるようにしてやると積分が0になってしまう問題が緩和されている。
何となくどんな重点関数を選べばよいのかわかってきたけれど、最適な重点関数の決め方はあるのだろうか?3.3.3節の説明によると理論的に最適な重点関数を決める定理があるらしいのだが、実用性は低いので上のように推定量の分散で収束性を判断するのが一般的のようだ。
次回は、Rによるモンテカルロ法入門に出てくる例と練習問題を解いていきたい。あともっとこまめに更新したい・・・