model
{
## first, standardise x's and coefficients##
for (j in 1 : p) {
   b[j] <- beta[j] / sd(x[ , j ])
   for (i in 1 : N) {z[i,j] <- (x[i,j] -  mean(x[,j])) / sd(x[,j])}
   }
   b0 <- beta0 - b[1] * mean(x[, 1]) - b[2] * mean(x[, 2]) - b[3] * mean(x[, 3])

# statistical model
for (i in 1 : N) {
   Y[i] ~ dnorm(mu[i], tau)
   mu[i] <- beta0 +beta[1] * z[i, 1] +  beta[2] * z[i, 2] +beta[3] * z[i, 3]
   }

# model validation
#  We only keep the J=ith element of each stat,
#  where ith datapoint is set to be NA for "leave-one-out" validation.

sresid_num <- (Ytrue[J] - mu[J]) 
Ey2<- sigma*sigma + mu[J]*mu[J]  
mui<-mu[J]
CPO <-sqrt(tau)* exp(-tau/2* (Ytrue[J] - mu[J])* (Ytrue[J] -mu[J]))
 
# priors
beta0 ~  dnorm(0, 0.00001)
for (j in 1 : p) {beta[j] ~ dnorm(0, 0.00001)}
sigma ~ dunif(0.01,100)         #  vague Gelman prior for sigma
tau<- 1/(sigma*sigma)

}  #  end of "exact" stacks model