###
# Create example data files for 20 independent chains
# with serial correlation of 0.25
###
set.seed(42)
tmpdir <- tempdir()
nsamples <- 1000
for(i in 1:20){
x <- matrix(nrow = nsamples+1, ncol=4)
colnames(x) <- c("alpha","beta","gamma", "nu")
x[,"alpha"] <- rnorm (nsamples+1, mean=0.025, sd=0.0025)^2
x[,"beta"] <- rnorm (nsamples+1, mean=53, sd=12)
x[,"gamma"] <- rbinom(nsamples+1, 20, p=0.25) + 1
x[,"nu"] <- rnorm (nsamples+1, mean=x[,"alpha"] * x[,"beta"], sd=1/x[,"gamma"])
# induce serial correlation of 0.25
x <- 0.75 * x[2:(nsamples+1),] + 0.25 * x[1:nsamples,]
write.table(
x,
file = file.path(
tmpdir,
paste("mcmc", i, "csv", sep=".")
),
sep = ",",
row.names = FALSE
)
}
# Read them back in as an mcmc.list object
data <- read.mcmc(
20,
file.path(tmpdir, "mcmc.#.csv"),
sep=",",
col.names=c("alpha","beta","gamma", "nu")
)
# Summary statistics
summary(data)
# Trace and Density Plots
plot(data)
# And check the necessary run length
mcgibbsit(data)
Run the code above in your browser using DataLab