水産資源学におけるベイズ統計の応用ワークショップのスライドの22,23枚目
- g(x)からM個のサンプルx1*, …, xM*を発生
- w*= f(x*)/g(x*)を計算
- x1*, …, xM*からw1*, …, wM*に比例する確率でサンプルをm個とってくる
- 得られた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))