setwd("K:/private/7440/Hatfield_mcmc_examples")
library(R2WinBUGS)
library(rbugs)
library(rjags)
library(BRugs)
library(mcmc)

## Generate toy data for analysis in multiple MCMC software options
## Simulated data and constants
N <- 100
set.seed(221983)
x <- rnorm(N,0,1)
beta <- c(4,2)
y <- rnorm(N,beta[1]+beta[2]*x,sqrt(2))
bprec <- c(1/10,1/10)
mu <- c(5,1)
prec <- 1/2

## Output data for fitting with WinBUGS, OpenBUGS, and JAGS directly
bugs.data(list(y=y,x=x,mu=mu,prec=prec,bprec=bprec,N=N),data.file="WB_data.txt")
dump(c('y','x','mu','prec','bprec','N'),file="JAGS_data.txt")

## Analyze using R2WinBUGS
r2fit <- bugs(
  data=list(y=y,x=x,mu=mu,prec=prec,bprec=bprec,N=N),
  inits=list(list(beta=c(0,0)),list(beta=c(1,1))),
  parameters.to.save='beta',model.file="WB_model.txt",
  bugs.directory="c:/Program Files (x86)/WinBUGS14/",
  n.chains=2,n.iter=10000,n.burnin=0,program="WinBUGS")
plot(r2fit)

## Analyze using rjags
jagsmodel <- jags.model(file="JAGS_model.txt",
  data=list(y=y,x=x,mu=mu,prec=prec,bprec=bprec,N=N),
 inits=list(list(beta=c(0,0)),list(beta=c(1,1))),n.chains=2)
jagsfit <- jags.samples(jagsmodel, variable.names='beta',n.iter=10000)
summary(jagsfit$beta,quantile,prob=c(.025,.5,.975))
par(mfrow=c(2,2))
plot(jagsfit$beta[1,,1],type='l',col='red',ylab='',main=expression(beta[0]))
lines(jagsfit$beta[1,,2],type='l',col='blue')
plot(jagsfit$beta[2,,1],type='l',col='red',ylab='',main=expression(beta[1]))
lines(jagsfit$beta[2,,2],type='l',col='blue')
plot(density(jagsfit$beta[1,,]),type='l',main='')
plot(density(jagsfit$beta[2,,]),type='l',main='')

## Analyze using BRugs
brugsfit <- BRugsFit(modelFile="WB_model.txt",
  data=list(y=y,x=x,mu=mu,prec=prec,bprec=bprec,N=N),
  inits=list(list(beta=c(0,0)),list(beta=c(1,1))),
  nIter=10000,nBurnin=0,
  numChains=2,parametersToSave='beta')
brugsfit$Stats
samplesHistory('beta')
samplesDensity('beta')

## Analyze using rbugs
rbugsfit <- rbugs(data=list(y=y,x=x,mu=mu,prec=prec,bprec=bprec,N=N),
  inits=list(list(beta=c(0,0)),list(beta=c(1,1))),
  paramSet='beta', model='OB_model.txt',
  n.chains=2,n.iter=10000,n.burnin=0,n.thin=1,
  bugs="C:/Program Files (x86)/OpenBUGS/OpenBUGS321/OpenBUGS.exe",
  OpenBugs=T,bugsWorkingDir=getwd())
par(mfrow=c(2,2))
plot(rbugsfit$chain1$'beta[1]',type='l',col='red',main=expression(beta[0]))
lines(rbugsfit$chain2$'beta[1]',col='blue')
plot(rbugsfit$chain1$'beta[2]',type='l',col='red',main=expression(beta[1]))
lines(rbugsfit$chain2$'beta[2]',col='blue')
plot(density(c(rbugsfit$chain1$'beta[1]',rbugsfit$chain2$'beta[1]')),type='l',main='')
plot(density(c(rbugsfit$chain1$'beta[2]',rbugsfit$chain2$'beta[2]')),type='l',main='')

## Analyze using mcmc
myll <- function(params,data){
  meen <- params[1] + params[2]*data$x
  ll <- sum(dnorm(data$y,meen,data$sigy,log=T))+sum(dnorm(params,data$mu,data$sigb,log=T))
 return(ll)
}
metropfit <- metrop(myll,initial=c(0,0),nbatch=10000,blen=1,scale=.2,
  data=list(x=x,y=y,mu=mu,sigy=sqrt(1/prec),sigb=sqrt(1/bprec)))
metropfit2 <- metrop(myll,initial=c(0,0),nbatch=10000,blen=1,scale=.2,
  data=list(x=x,y=y,mu=mu,sigy=sqrt(1/prec),sigb=sqrt(1/bprec)))
metropfit$accept
metropfit2$accept
first <- 101
apply(rbind(metropfit$batch[first:10000,],metropfit2$batch[first:10000,]),
  2,quantile,prob=c(.025,.5,.975))
par(mfrow=c(2,2))
plot(metropfit$batch[first:10000,1],type='l',col='red',main=expression(beta[0]))
lines(metropfit2$batch[first:10000,1],col='blue')
plot(metropfit$batch[first:10000,2],type='l',col='red',main=expression(beta[1]))
lines(metropfit2$batch[first:10000,2],col='blue')
plot(density(c(metropfit$batch[first:10000,1],metropfit2$batch[first:10000,1])),type='l',main='')
plot(density(c(metropfit$batch[first:10000,2],metropfit2$batch[first:10000,2])),type='l',main='')