library(R2WinBUGS)
library(MASS)  # need to mvrnorm
library(MCMCpack) # need for rwish

# Generate synthetic data
N = 500
S = matrix(c(1,.2,.2,5),nrow=2)
y = mvrnorm(N,c(1,3),S)


# Set up for WinBUGS
mu0 = as.vector(c(0,0))
S2 = matrix(c(1,0,0,1),nrow=2)/1000
S3 =  matrix(c(1,0,0,1),nrow=2)/10000
data=list("y","N","S2","S3","mu0")
inits=function(){list( mu=mvrnorm(1,mu0,matrix(c(10,0,0,10),nrow=2) ), 
tau = rwish(3,matrix(c(.02,0,0,.04),nrow=2)) )}

# Run WinBUGS
multi_norm.sim = bugs(data,inits,model.file="mult_normal.bug",
parameters=c("mu","tau"),n.chains = 2,n.iter=4010,n.burnin=10,n.thin=1,
bugs.directory="c:/Program Files/WinBUGS14/",codaPkg=FALSE)

print(multi_norm.sim,digits=3)
print(solve(S)) # prints the S^{-1} = tau

windows()
plot(multi_norm.sim)

windows()
par(mfrow=c(3,2))
for (i in 1:6)
{
ts.plot(multi_norm.sim$sims.array[,1,i])
}

windows()
par(mfrow=c(3,2))
for (i in 1:6)
{
acf(multi_norm.sim$sims.array[,1,i])
}