通过 R 在 WinBUGS 中设置.seed答案

【问题标题】：set.seed in WinBUGS through R通过 R 在 WinBUGS 中设置.seed
【发布时间】：2015-06-13 15:51:21
【问题描述】：

是否可以在WinBUGS 中设置种子以重现参数估计，就像在R 中使用set.seed 一样？

【问题讨论】：

【解决方案1】：

下面的代码建议您可以通过在每次运行 WinBUGS 模型之前立即设置种子来重现 WinBUGS 到 R 中的估计值。

前四个模型运行之前紧跟相同的set.seed 语句。最后两个模型运行不是。根据all.equal 声明，前四个模型运行返回相同的估计值。最后两个模型运行没有。

####################################################################################
####################################################################################


library(R2WinBUGS)

n <- 15
x <- 1:15
y <- c(32.46, 38.38, 40.92, 22.27, 34.64, 33.53, 26.62, 25.26, 23.67, 20.54, 21.11, 17.00, 16.61, 19.32, 22.29)


print(summary(lm(y ~ x)))


####################################################################################
####################################################################################


set.seed(1234)


sink("linreg.txt")
cat("
model {

# Priors
 alpha ~ dnorm(0,0.001)
 beta  ~ dnorm(0,0.001)
 sigma ~ dunif(0,  100)
 tau <- 1/ (sigma * sigma)

# Likelihood
 for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau) 
    mu[i] <- alpha + beta*x[i]
 }

}
",fill=TRUE)
sink()


win.data <- list("x","y", "n")

inits <- function(){list(alpha=rnorm(1), beta=rnorm(1), sigma = rlnorm(1))}

params <- c("alpha", "beta", "sigma")

nc =    2  
ni = 1000  
nb =  500
nt =    5

out1 <- bugs(data = win.data, inits = inits, parameters = params, 
            model = "linreg.txt", bugs.directory="c:/WinBUGS14/",
            n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni)

print(out1, dig = 4)


####################################################################################
####################################################################################


set.seed(1234)


sink("linreg.txt")
cat("
model {

# Priors
 alpha ~ dnorm(0,0.001)
 beta  ~ dnorm(0,0.001)
 sigma ~ dunif(0,  100)
 tau <- 1/ (sigma * sigma)

# Likelihood
 for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau) 
    mu[i] <- alpha + beta*x[i]
 }

}
",fill=TRUE)
sink()


win.data <- list("x","y", "n")

inits <- function(){list(alpha=rnorm(1), beta=rnorm(1), sigma = rlnorm(1))}

params <- c("alpha", "beta", "sigma")

nc =    2  
ni = 1000  
nb =  500
nt =    5

out2 <- bugs(data = win.data, inits = inits, parameters = params, 
            model = "linreg.txt", bugs.directory="c:/WinBUGS14/",
            n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni)

print(out2, dig = 4)


####################################################################################
####################################################################################


set.seed(1234)


sink("linreg.txt")
cat("
model {

# Priors
 alpha ~ dnorm(0,0.001)
 beta  ~ dnorm(0,0.001)
 sigma ~ dunif(0,  100)
 tau <- 1/ (sigma * sigma)

# Likelihood
 for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau) 
    mu[i] <- alpha + beta*x[i]
 }

}
",fill=TRUE)
sink()


win.data <- list("x","y", "n")

inits <- function(){list(alpha=rnorm(1), beta=rnorm(1), sigma = rlnorm(1))}

params <- c("alpha", "beta", "sigma")

nc =    2  
ni = 1000  
nb =  500
nt =    5

out3 <- bugs(data = win.data, inits = inits, parameters = params, 
            model = "linreg.txt", bugs.directory="c:/WinBUGS14/",
            n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni)

print(out3, dig = 4)


####################################################################################
####################################################################################


set.seed(1234)


sink("linreg.txt")
cat("
model {

# Priors
 alpha ~ dnorm(0,0.001)
 beta  ~ dnorm(0,0.001)
 sigma ~ dunif(0,  100)
 tau <- 1/ (sigma * sigma)

# Likelihood
 for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau) 
    mu[i] <- alpha + beta*x[i]
 }

}
",fill=TRUE)
sink()


win.data <- list("x","y", "n")

inits <- function(){list(alpha=rnorm(1), beta=rnorm(1), sigma = rlnorm(1))}

params <- c("alpha", "beta", "sigma")

nc =    2  
ni = 1000  
nb =  500
nt =    5

out4 <- bugs(data = win.data, inits = inits, parameters = params, 
            model = "linreg.txt", bugs.directory="c:/WinBUGS14/",
            n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni)

print(out4, dig = 4)


####################################################################################
####################################################################################


sink("linreg.txt")
cat("
model {

# Priors
 alpha ~ dnorm(0,0.001)
 beta  ~ dnorm(0,0.001)
 sigma ~ dunif(0,  100)
 tau <- 1/ (sigma * sigma)

# Likelihood
 for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau) 
    mu[i] <- alpha + beta*x[i]
 }

}
",fill=TRUE)
sink()


win.data <- list("x","y", "n")

inits <- function(){list(alpha=rnorm(1), beta=rnorm(1), sigma = rlnorm(1))}

params <- c("alpha", "beta", "sigma")

nc =    2  
ni = 1000  
nb =  500
nt =    5

out5 <- bugs(data = win.data, inits = inits, parameters = params, 
            model = "linreg.txt", bugs.directory="c:/WinBUGS14/",
            n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni)

print(out5, dig = 4)


####################################################################################
####################################################################################


sink("linreg.txt")
cat("
model {

# Priors
 alpha ~ dnorm(0,0.001)
 beta  ~ dnorm(0,0.001)
 sigma ~ dunif(0,  100)
 tau <- 1/ (sigma * sigma)

# Likelihood
 for (i in 1:n) {
    y[i] ~ dnorm(mu[i], tau) 
    mu[i] <- alpha + beta*x[i]
 }

}
",fill=TRUE)
sink()


win.data <- list("x","y", "n")

inits <- function(){list(alpha=rnorm(1), beta=rnorm(1), sigma = rlnorm(1))}

params <- c("alpha", "beta", "sigma")

nc =    2  
ni = 1000  
nb =  500
nt =    5

out6 <- bugs(data = win.data, inits = inits, parameters = params, 
            model = "linreg.txt", bugs.directory="c:/WinBUGS14/",
            n.thin = nt, n.chains = nc, n.burnin = nb, n.iter = ni)

print(out6, dig = 4)

####################################################################################
####################################################################################

all.equal(out1, out2)
all.equal(out1, out3)
all.equal(out1, out4)
all.equal(out1, out5)
all.equal(out1, out6)

【讨论】：