インポータンス・サンプリングの基礎 - ryamadaの遺伝学・遺伝統計学メモ

あるドメインに生起確率分布があるとする
そのドメインに、値があるとして、生起確率を考慮した期待値の推定をする
このとき、モンテカルロで推定するとする
どこもかしこも均一にランダム発生して値を計算してその平均をとるより
値確認をする頻度にバラツキを入れる方が、精度がよくなる場合があること
そのようなバラツキを入れたサンプリング(重みづけサンプリング)をするにあたっては、「もともとの生起確率分布」に応じてサンプリングするとよい
問題は、そのような分布がわかっていない場合もあることだが、それでも、何かしらとっかかりがあるなら、重みづけサンプリングがよい
以下の例では次のように設定で、このことをシミュレーションしている
xが[0,1]の範囲をとる
ベータ分布p(x)に従ってxは分布し
それぞれのxに対して、 $f(x)=x^2+3$ のような値をとるのだが、じゃあ、分布p(x)の下でのf(x)の期待値は？という問題
初めに、x:[0,1]に一様乱数を発生させ、 $p(x)f(x)$ の期待値を出す
次に、p(x)はわからないのだけれど、適当にx乱数を発生させることにする。ここでは「パラメタ値を適当に変えてベータ乱数」を発生させることにする。この乱数の分布をq(x)とする
そのようにして得られたx乱数に関して、f(x)を計算し、重みとして $p(x)/q(x)$ を与えて、重みつき平均を出す
こうすると、期待値は、どれでやっても、それなりだが、分散はp(x)=q(x)のときに小さくなることが読み取れる
ここで、乱数発生分布にベータ分布を使ったけれど、これを適当に変えて、もとの分布p(x)と違っていてもよい(その方がふつうか…)

a <- 2
b <- 3
my.f <- function(x){
	x^2+3
	#matrix(1,nrow(x),ncol(x))
}
my.p <- function(x,a=2,b=3){
	dbeta(x,a,b)
}

m <- 10^2

n.iter <- 10^3

x <- matrix(runif(m*n.iter),m,n.iter)
fx <- my.f(x)
px <- my.p(x,a,b)

fpx <- fx*px
mu <- apply(fpx,2,mean)
hist(mu)


As <- seq(from=0.5,to=5,length=20)
Bs <- seq(from=0.5,to=5,length=20)

ABs <- expand.grid(As,Bs)
mus <- matrix(0,length(ABs[,1]),n.iter)
for(i in 1:length(ABs[,1])){
	x <- matrix(rbeta(m*n.iter,ABs[i,1],ABs[i,2]),m,n.iter)
	px <- my.p(x,a,b)
	qx <- my.p(x,ABs[i,1],ABs[i,2])
	wx <- px/qx
	fx <- my.f(x)
	fwx <- fx*wx
	#mu. <- apply(fwx,2,sum)/apply(wx,2,sum)
	mu. <- apply(fwx,2,mean)
	mus[i,] <- mu.
}
library(rgl)
plot3d(ABs[,1],ABs[,2],log(apply(mus,1,var)))