COV
切り捨てられたベータ分布と一致する の値をシミュレートしたいだけの場合はCOV
、データのリストから省略できます。例えば:
library(R2jags)
cat('
model {
for (i in 1:n.sites) {
COV[i] ~ dbeta(a[i], b[i])T(r1[i], r2[i])
}
}', file={M <- tempfile()})
dat <- list(n.sites=5, a=c(7.1, 2.2, 13, 10, 13), b=c(25, 11, 44, 27, 44),
r1=c(0.05, 0.1, 0.2, 0.1, 0.2), r2=c(0.15, 0.3, 0.6, 0.3, 0.6))
j <- jags(dat, NULL, 'COV', M, 1, 10000, DIC=FALSE, n.burnin=0, n.thin=1)
これは、10000 個のシミュレートされたベクトルのマトリックスの上位数行COV
です...
head(j$BUGSoutput$sims.list$COV)
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.1169165 0.2889155 0.2543063 0.2083161 0.2426788
## [2,] 0.1494647 0.1430956 0.2867575 0.2410594 0.2795923
## [3,] 0.1200414 0.2093230 0.2736719 0.2189734 0.2469634
## [4,] 0.1472082 0.1442609 0.2911482 0.2625216 0.2714883
## [5,] 0.1403574 0.1100977 0.2556352 0.1918480 0.2353231
## [6,] 0.1310404 0.1677148 0.3011752 0.1974136 0.2131811
編集
既知の分布からサンプリングしているだけなので、純粋な R でこれを行うこともできます (distr
パッケージは、切り捨てられた分布からのサンプリングに役立ついくつかの関数を提供します)。
library(distr)
n <- 10000 # How many samples?
COV <- mapply(function(shape1, shape2, min, max) {
r(Truncate(Beta(shape1, shape2), min, max))(n)
}, shape1=c(7.1, 2.2, 13, 10, 13), shape2=c(25, 11, 44, 27, 44),
min=c(0.05, 0.1, 0.2, 0.1, 0.2), max=c(0.15, 0.3, 0.6, 0.3, 0.6))
shape1
上記では、 、shape2
、min
およびの等しい長さのベクトルを渡す必要があります。これにより、それぞれに対してランダムなベータ分布変量が生成さmax
れ、対応する、 、およびが順番に生成されます。mapply
n
shape1
shape2
min
max
COV
(純粋な R サンプル)の列のカーネル密度をj$BUGSoutput$sims.list$COV
(JAGS サンプル)の列のカーネル密度と比較できます。
par(mfrow=c(3, 2), mar=c(3, 0.5, 0.5, 0.5))
sapply(1:5, function(i) {
djags <- density(j$BUGSoutput$sims.list$COV[, i])
dr <- density(COV[, i])
plot(djags, lwd=4, col='gray80', main='', ylab='', xlab='', yaxt='n',
ylim=c(0, max(djags$y, dr$y)))
lines(dr)
})
plot.new()
legend('topleft', c('JAGS', 'R'), col=c('gray80', 'black'), lwd=c(4, 1), bty='n')