运行 R2WinBUGS 中 n 个数据集的模拟研究
Run simulation studies in R2WinBUGS for n datasets
我在弄清楚如何使用 R2WinBUGS 运行 进行一些模拟研究时遇到问题。目的是模拟 n 个数据集(目标是 1000,但从 10 个开始),并将它们作为矩阵全部放入 R2WinBUGS 代码中,以便当它移植到 WinBUGS 时,它将 运行 产生估计n 个数据集。这是我目前拥有的:
模特:
model{
alpha0 ~ dnorm(66.6, 0.01)
alpha1 ~ dnorm(0.3, 0.01)
alpha2 ~ dnorm(100, 0.01)
alpha3 ~ dnorm(0.2, 0.01)
beta0 ~ dnorm(35, 0.01)
beta1 ~ dnorm(80, 0.01)
tau ~ dgamma(0.3,1)
for(k in 1:Ndat) {
y[k] ~ dnorm(mu[k], tau)
mu[k] <- ((alpha0/(1 + exp(-alpha1*(28-beta0)))) + (alpha2/(1 + exp(-alpha3*(28-beta1)))))
}
}
我使用的bug代码是:
grapedat.sim = bugs(data = list('Ndat' = Ndat, 'y' = p.y[,1]),inits,
model.file="H:/R coding/R2WinBUGS/multsimt1.bug",
parameters=c("alpha0","alpha1","alpha2","alpha3","beta0","beta1","tau"),
n.chains=1,n.iter=8000,n.sim = 6000,
n.burnin=2000,n.thin=1,
bugs.directory="H:/WinBUGS14",
codaPkg=FALSE,
debug = T)
其中 Ndat 是数据集的数量,p.y。是一个13 x n的矩阵,inits是:
inits <- function(){
list(alpha0=70, alpha1=0.4, tau=0.15, alpha2=105, alpha3=0.25,beta0 =
40, beta1 = 85)
}
有什么帮助吗?
理论上,您可以在单个 BUGS 模型中执行此操作,方法是使用第二个(外部)循环,然后将 y 和 mu 索引为矩阵,并将系数作为向量,例如(未经测试):
model{
for(d in 1:n){
alpha0[d] ~ dnorm(66.6, 0.01)
alpha1[d] ~ dnorm(0.3, 0.01)
alpha2[d] ~ dnorm(100, 0.01)
alpha3[d] ~ dnorm(0.2, 0.01)
beta0[d] ~ dnorm(35, 0.01)
beta1[d] ~ dnorm(80, 0.01)
tau[d] ~ dgamma(0.3,1)
for(k in 1:Ndat) {
y[k,d] ~ dnorm(mu[k,d], tau[d])
mu[k,d] <- ((alpha0[d]/(1 + exp(-alpha1[d]*(28-beta0[d])))) +
(alpha2[d]/(1 + exp(-alpha3[d]*(28-beta1[d])))))
}
}
}
这对于 n=10 是可行的,但对于 n=1000 将非常笨拙,因此通常不是解决问题的好方法。相反,我会 运行 每个 dataset/model 独立,在 R 中循环。如果您愿意使用 JAGS 而不是 WinBUGS,那么您可以在 运行jags 包 - 即阅读本文第 2.9 节:https://www.jstatsoft.org/article/view/v071i09
马特
我在弄清楚如何使用 R2WinBUGS 运行 进行一些模拟研究时遇到问题。目的是模拟 n 个数据集(目标是 1000,但从 10 个开始),并将它们作为矩阵全部放入 R2WinBUGS 代码中,以便当它移植到 WinBUGS 时,它将 运行 产生估计n 个数据集。这是我目前拥有的:
模特:
model{
alpha0 ~ dnorm(66.6, 0.01)
alpha1 ~ dnorm(0.3, 0.01)
alpha2 ~ dnorm(100, 0.01)
alpha3 ~ dnorm(0.2, 0.01)
beta0 ~ dnorm(35, 0.01)
beta1 ~ dnorm(80, 0.01)
tau ~ dgamma(0.3,1)
for(k in 1:Ndat) {
y[k] ~ dnorm(mu[k], tau)
mu[k] <- ((alpha0/(1 + exp(-alpha1*(28-beta0)))) + (alpha2/(1 + exp(-alpha3*(28-beta1)))))
}
}
我使用的bug代码是:
grapedat.sim = bugs(data = list('Ndat' = Ndat, 'y' = p.y[,1]),inits,
model.file="H:/R coding/R2WinBUGS/multsimt1.bug",
parameters=c("alpha0","alpha1","alpha2","alpha3","beta0","beta1","tau"),
n.chains=1,n.iter=8000,n.sim = 6000,
n.burnin=2000,n.thin=1,
bugs.directory="H:/WinBUGS14",
codaPkg=FALSE,
debug = T)
其中 Ndat 是数据集的数量,p.y。是一个13 x n的矩阵,inits是:
inits <- function(){
list(alpha0=70, alpha1=0.4, tau=0.15, alpha2=105, alpha3=0.25,beta0 =
40, beta1 = 85)
}
有什么帮助吗?
理论上,您可以在单个 BUGS 模型中执行此操作,方法是使用第二个(外部)循环,然后将 y 和 mu 索引为矩阵,并将系数作为向量,例如(未经测试):
model{
for(d in 1:n){
alpha0[d] ~ dnorm(66.6, 0.01)
alpha1[d] ~ dnorm(0.3, 0.01)
alpha2[d] ~ dnorm(100, 0.01)
alpha3[d] ~ dnorm(0.2, 0.01)
beta0[d] ~ dnorm(35, 0.01)
beta1[d] ~ dnorm(80, 0.01)
tau[d] ~ dgamma(0.3,1)
for(k in 1:Ndat) {
y[k,d] ~ dnorm(mu[k,d], tau[d])
mu[k,d] <- ((alpha0[d]/(1 + exp(-alpha1[d]*(28-beta0[d])))) +
(alpha2[d]/(1 + exp(-alpha3[d]*(28-beta1[d])))))
}
}
}
这对于 n=10 是可行的,但对于 n=1000 将非常笨拙,因此通常不是解决问题的好方法。相反,我会 运行 每个 dataset/model 独立,在 R 中循环。如果您愿意使用 JAGS 而不是 WinBUGS,那么您可以在 运行jags 包 - 即阅读本文第 2.9 节:https://www.jstatsoft.org/article/view/v071i09
马特