Is it possible to set a seed in WinBUGS to reproduce parameter estimates like one can do with set.seed in R?
Asked
Active
Viewed 632 times
2
Mark Miller
- 12,483
- 23
- 78
- 132
1 Answers
3
The code below suggests that you can reproduce estimates in WinBUGS through R by setting a seed immediately before each run of the WinBUGS model.
The first four model runs are immediately preceded by the same set.seed statement. The last two model runs are not. According to the all.equal statement the first four model runs return the identical estimates. The last two model runs do not.
####################################################################################
####################################################################################
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)
Mark Miller
- 12,483
- 23
- 78
- 132