# Code written by Eric Ward, eric.ward@noaa.gov on 07/01/08
# This script uses either OpenBUGS or WinBUGS to do estimation of 
# the multiple linear regression model with constant variance.
library("R2WinBUGS")

# covariates consist of predictors X (n x p, including 1s)
# response is univariate Y (n x 1)
# 
n = 10
p = 3
# Create simulated data for the example
X <- matrix(cbind(rep(1,n),seq(1,n),runif(n)),nrow=n,ncol=p)
Y <- rnorm(n= n, mean = X%*%c(0.7,0.3,0.1),1)
# Assign the prior means/variances to the 
B.prior_means <- rep(0,p)
B.prior_vars <- rep(10,p)
# Assign the shape / scale parameters
# Note: I'm using the ScaledInvChi2 distribution, following (Gelman 2004).
# Shape (v) interpreted as the amount of information the prior contributes,
# Scale (theta) interpreted as prior variance
# sig2 ~ ScaledInvChi2 (v, theta) equivalent to sig2 ~ InvGamma(v/2,theta/2)
# or tau = (1/sig2) ~ Gamma(shape = v/2, rate= v/2s2) with E[X] = 1/s2
sigma.v = 1
sigma.theta = 1
tau.prior = c(sigma.v/2, (sigma.v/2)*sigma.theta)

# this is the model string, generalized for ANY dataset
modelString = ("
model {
  for(i in 1:p) {
     B.tau[i] <- 1/B.prior_vars[i]
     B[i] ~ dnorm(B.prior_means[i],B.tau[i]);
  }
  tau ~ dgamma(tau.prior[1],tau.prior[2]);
  sigma <- sqrt(1/tau);
  for(i in 1:n) {
    pred[i] <- inprod(X[i,1:p],B[1:p]);
    Y[i] ~ dnorm(pred[i],tau);
  } 
}
")

burn = 10000
chainLength = 10000
write(modelString,file="model.txt",append=FALSE)
data <- list("Y","X","n","p","B.prior_means","B.prior_vars","tau.prior");
parameters <- c("B","sigma")
mod <- NA
try(mod <- bugs (data,inits=NULL,parameters,DIC=TRUE,"model.txt",n.chains=1,n.burnin=burn, 
   n.thin=1,n.iter=(burn+chainLength),program="openbugs",debug=FALSE),silent=TRUE)
attach.bugs(mod,overwrite=TRUE)

# look at model summary
mod
# plot parameters
par(mfrow=c(2,2))
for(i in 1:p) hist(B[,i],100,xlab=paste("B",i-1,sep=""),main="")
hist(sigma,100,xlab=expression(sigma),main="")