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
  for (i in 1:N) {
      Y[i] ~ ddexp(mu[i],tau);       # DE errors
      mu[i] <- beta0 + beta[1]*z[i,1]+beta[2]*z[i,2]+beta[3]*z[i,3];
      ds[i] <- 2*tau*(max(Y[i],mu[i])-min(Y[i],mu[i]))
               - 2*log(tau/2);                             # DE errors
  }
# 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);
}