model
{
# 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]) ;
      }
  }
# Model
  PI <- 3.141593;
  d <- 4;  # degrees of freedom for t
  d2 <- d/2;  # half of that
  for (i in 1:N) {
      prec[i] <- w[i]*tau;           # t_d errors via normal scale mixture
      Y[i] ~ dnorm(mu[i],prec[i]);
      w[i] ~ dgamma(d2,d2);
      mu[i] <- beta0 + beta[1]*z[i,1]+beta[2]*z[i,2]+beta[3]*z[i,3];
      ds[i] <- w[i]*tau*pow((Y[i]-mu[i]),2) 
               - log((w[i]*tau)/(2*PI));  # t_d errors via normal scale mixture
  }
# Priors 
  beta0 ~  dnorm(0,.00001);
  for (j in 1:p) {
     beta[j] ~ dnorm(0,.00001); # coeffs independent
  }
  tau ~ dgamma(1.0E-3,1.0E-3);
}