RUNJAGS-在没有事先进行模拟的情况下设置种子
RUNJAGS-set seed without prior for simulations
我正在使用 runjags 从正态分布中抽取一些数据。对于我用于模拟的参数,我没有任何先验知识。 runjages 似乎没有使用参数来固定种子:list(".RNG.name"="base::Super-Duper", ".RNG.seed"=1)
。
我将参数更改为 list(muOfClustsim=rep(1, npop), ".RNG.name"="base::Super-Duper", ".RNG.seed"=1)
但它也不起作用。有没有办法在 runjags 中修复此类模型的种子?
这是一个最小的可重现示例:
library(runjags)
npop=3
nrep=10
sdpop=7
sigma=5
seed=4
set.seed(seed)
N = npop*nrep # nb of observations
## Population identity of each individual used to sample genotypes but not used for common garden test
pop <- rep(1:npop, each=nrep)
muOfClustsim <- rnorm(npop, 0, sdpop) # vector of population means
(tausim <- 1/(sigma*sigma)) # precision of random individual error
# parameters are treated as data for the simulation step
data <- list(N=N, pop=pop, muOfClustsim=muOfClustsim, tausim=tausim)
## JAG model
txtstring <- "
data{
# Likelihood:
for (i in 1:N){
ysim[i] ~ dnorm(eta[i], tausim) # tau is precision (1 / variance)
eta[i] <- muOfClustsim[pop[i]]
}
}
model{
fake <- 0
}
"
## Initial values with seed for reproducibility
initssim <- list(".RNG.name"="base::Super-Duper", ".RNG.seed"=1)
##initssim <- list(muOfClustsim=rep(1, npop), ".RNG.name"="base::Super-Duper", ".RNG.seed"=1)
## Simulate with jags
set.seed(seed)
out <- run.jags(txtstring, data = data, monitor=c("ysim"), sample=1, n.chains=1, inits=initssim, summarise=FALSE)
## reformat the outputs
(ysim1 <- coda::as.mcmc(out)[1:N])
set.seed(seed)
out <- run.jags(txtstring, data = data, monitor=c("ysim"), sample=1, n.chains=1, inits=initssim, summarise=FALSE)
## reformat the outputs
(ysim2 <- coda::as.mcmc(out)[1:N])
identical(ysim1, ysim2)
.RNG.name / .RNG.seed / .RNG.state 初始值适用于模型(或更具体地说,模型中的链),并且不在模型中使用数据块。这意味着(据我所知)没有办法在 JAGS <= 4.3 中重现数据块内的任何随机表示。这是可以为 JAGS 的未来版本添加的东西,但我担心它的优先级较低,因为在将数据传递给 JAGS 之前总是可以(通常更好)在 R 中模拟数据。
在你的情况下,答案(假设你想使用 JAGS)是在模型块而不是数据块中进行模拟:
txtstring <- "
model{
# Likelihood:
for (i in 1:N){
ysim[i] ~ dnorm(eta[i], tausim) # tau is precision (1 / variance)
eta[i] <- muOfClustsim[pop[i]]
}
}
"
然后您的其余代码将按预期运行 §。值得注意的是,数据模拟通常是一项比 JAGS 更适合 R 的任务,但我假设在这种情况下您有特定的原因想要使用 JAGS...
马特
§ 尽管您通常不应期望双打严格相等,例如:
identical(0.1+0.2, 0.3)
但是:
abs(0.3 - (0.1+0.2)) < sqrt(.Machine$double.eps)
甚至更好:
isTRUE(all.equal(0.1+0.2, 0.3))
我正在使用 runjags 从正态分布中抽取一些数据。对于我用于模拟的参数,我没有任何先验知识。 runjages 似乎没有使用参数来固定种子:list(".RNG.name"="base::Super-Duper", ".RNG.seed"=1)
。
我将参数更改为 list(muOfClustsim=rep(1, npop), ".RNG.name"="base::Super-Duper", ".RNG.seed"=1)
但它也不起作用。有没有办法在 runjags 中修复此类模型的种子?
这是一个最小的可重现示例:
library(runjags)
npop=3
nrep=10
sdpop=7
sigma=5
seed=4
set.seed(seed)
N = npop*nrep # nb of observations
## Population identity of each individual used to sample genotypes but not used for common garden test
pop <- rep(1:npop, each=nrep)
muOfClustsim <- rnorm(npop, 0, sdpop) # vector of population means
(tausim <- 1/(sigma*sigma)) # precision of random individual error
# parameters are treated as data for the simulation step
data <- list(N=N, pop=pop, muOfClustsim=muOfClustsim, tausim=tausim)
## JAG model
txtstring <- "
data{
# Likelihood:
for (i in 1:N){
ysim[i] ~ dnorm(eta[i], tausim) # tau is precision (1 / variance)
eta[i] <- muOfClustsim[pop[i]]
}
}
model{
fake <- 0
}
"
## Initial values with seed for reproducibility
initssim <- list(".RNG.name"="base::Super-Duper", ".RNG.seed"=1)
##initssim <- list(muOfClustsim=rep(1, npop), ".RNG.name"="base::Super-Duper", ".RNG.seed"=1)
## Simulate with jags
set.seed(seed)
out <- run.jags(txtstring, data = data, monitor=c("ysim"), sample=1, n.chains=1, inits=initssim, summarise=FALSE)
## reformat the outputs
(ysim1 <- coda::as.mcmc(out)[1:N])
set.seed(seed)
out <- run.jags(txtstring, data = data, monitor=c("ysim"), sample=1, n.chains=1, inits=initssim, summarise=FALSE)
## reformat the outputs
(ysim2 <- coda::as.mcmc(out)[1:N])
identical(ysim1, ysim2)
.RNG.name / .RNG.seed / .RNG.state 初始值适用于模型(或更具体地说,模型中的链),并且不在模型中使用数据块。这意味着(据我所知)没有办法在 JAGS <= 4.3 中重现数据块内的任何随机表示。这是可以为 JAGS 的未来版本添加的东西,但我担心它的优先级较低,因为在将数据传递给 JAGS 之前总是可以(通常更好)在 R 中模拟数据。
在你的情况下,答案(假设你想使用 JAGS)是在模型块而不是数据块中进行模拟:
txtstring <- "
model{
# Likelihood:
for (i in 1:N){
ysim[i] ~ dnorm(eta[i], tausim) # tau is precision (1 / variance)
eta[i] <- muOfClustsim[pop[i]]
}
}
"
然后您的其余代码将按预期运行 §。值得注意的是,数据模拟通常是一项比 JAGS 更适合 R 的任务,但我假设在这种情况下您有特定的原因想要使用 JAGS...
马特
§ 尽管您通常不应期望双打严格相等,例如:
identical(0.1+0.2, 0.3)
但是:
abs(0.3 - (0.1+0.2)) < sqrt(.Machine$double.eps)
甚至更好:
isTRUE(all.equal(0.1+0.2, 0.3))