drkcore

2009/08/18 22:13:35

SIR(sampling-importance-resampling)

水産資源学におけるベイズ統計の応用ワークショップのスライドの22,23枚目

1. g(x)からM個のサンプルx1*, …, xM*を発生
2. w*= f(x*)/g(x*)を計算
3. x1*, …, xM*からw1*, …, wM*に比例する確率でサンプルをm個とってくる
4. 得られたm個のサンプルはf(x)からの事後サンプル

となっているのだが、23枚目のRのコードは

n2 <-sample(length(x1), 1000, replace=T, prob=dlnorm(x1,1.1,0.6))

とwを求めずにいきなりサンプリングしているので？となったが、一様分布だから重みを計算する必要ないのね。

わざわざ書くならこうか。

x1 <-runif(100000,0,60)
w <- dlnorm(x1,1.1,0.6) / (dunif(x1,0,60))
n2 <-sample(length(x1), 1000, replace=T, prob=w)
x2 <-x1[n2]
plot(density(x2))

Comments