重点サンプリング (4)（2014/10/11）のつづき。今回は、Rによるモンテカルロ法入門のテキストの例3.5を解いてみたい。

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

の積分を重点サンプリングで計算しろというもの。似たような積分は (1)～(3) までですでに取り上げてきた。そこでは、標準正規分布からのサンプリングでは4.5以上のサンプルを得ることはほぼありえないため、重点関数として平均値をずらした正規分布を用いてサンプリングした。今回は、重点関数として平均値をずらした正規分布ではなく、4.5で切り詰めた指数分布を使えというのが課題になる。

切断指数分布

切断指数分布（truncated exponential distribution）という分布は初めて聴いたのだが、指数分布から4.5以上の区間を取り出してきて、その区間の面積が1になるように再調整した分布ということがわかった。厳密な定義は、WikipediaのTruncated distributionの項を参照のこと。指数分布を4.5で切り詰めた分布は、

$\displaystyle g(y) = \frac{e^{-y}}{\int_{4.5}^{\infty} e^{-x} dx} = e^{-(y-4.5)}$

で表される。Wikipediaの定義からもちゃんと導出できることは確認した。条件として4.5より小さい範囲では、 $g(y) = 0$ となるので注意。早速、どのような分布なのか確認してみよう。

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]

f:id:aidiary:20141011124535p:plain

全区間で積分すると1になることも確認できた。

切断指数分布を用いた重点サンプリング

次に重点関数として4.5で切り詰めた切断指数分布を用いてモンテカルロ積分を計算してみよう。ついでに重点サンプリング (2)（2014/9/21）で用いた正規分布の結果もいっしょに収束性の違いを評価してみた。青色が本来のサンプリング分布である標準正規分布。緑色が1つめの重点関数候補である平均をずらした正規分布。赤色が2つめの重点関数候補である切断指数分布。

f:id:aidiary:20141011124538p:plain

切断指数分布の関数は、引数に 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が正規分布から生成できないためほとんどのサンプルで $h(x) = 0$ となってしまい積分値が0になってしまう。重点サンプリングを使うとその問題がなくなる。

f:id:aidiary:20141011124540p:plain

重点サンプリングで、重点関数に正規分布を使った場合と切断指数分布を使った場合を比較すると、切断指数分布を使った方が収束性が早く、分散も小さいことがわかる。おそらく、正規分布では4.5より小さい左半分のサンプルが無駄になるのに対し、切断指数分布では必ず4.5以上のサンプルが生成され、無駄になるサンプルがないためだと考えられる。