Example of importance sampling

マルコフ連鎖モンテカルロ法 (統計ライブラリー)の1.2.3を参考にしつつ。

rejection samplingでは、target distributionからの乱数を間接的に取ってきているわけですが、今度はxの関数h(x)の期待値に興味があるとしましょう。で、積分が難しいので、サンプリングしてどうこうしまようというのがこのimportance sampling。ここではh(x)=xとして、普通に期待値を計算させてきましょう。g(x)はh(x)p(x)になるべく近く、サンプリング可能であるとする、ということなので、またガンマ関数から取ってこれないとして、コーシー分布から取ってくるというのをやります。定数Mで覆えるというとかいうのはここではありません。

plot(function(x){dgamma(x,shape=10)},0,30)
plot(function(x){dcauchy(x,scale=2.5,location=10)},0,30,add=TRUE)

f:id:syou6162:20161024213125p:plain

> t <- 1000
> x <- rcauchy(t,scale=2.5,location=10)
> w <- dgamma(x,shape=10) / (dcauchy(x,scale=2.5,location=10))
> sum(w*x) / t
[1] 10.20139