#############
#  Here is WinBUGS code for normal scale mixture model
#  (Carlin and Louis, 2008, Example 4.7)
#############

model
{
  for(i in 1:N) {
    # Duplicate the data
    Y1[i] <- Y[i]
    Y2[i] <- Y[i]
    Y3[i] <- Y[i]

    # Weighted precision parameters
    tau1[i] <- tau0[1]/lambda1[i];
    tau2[i] <- tau0[2]/lambda2[i];
    tau3[i] <- tau0[3]/lambda3[i];

    # Mean structures (all same)
    Y1[i] ~ dnorm(theta[1],tau1[i]);
    Y2[i] ~ dnorm(theta[2],tau2[i]);
    Y3[i] ~ dnorm(theta[3],tau3[i]);

    # Error distributions
    lambda1[i] <- 1;                # M1 = e_i ~ N(0,tau1)
    lambda2[i] <- 1/inv.lambda2[i]; # So that lambda ~ IG
    inv.lambda2[i] ~ dgamma(1,1);   # M2 = e_i ~ t_2(0, tau2)
    lambda3[i] ~ dexp(0.5);         # M3 = e_i ~ DE(0,tau3)
  }
  # Priors
  for (k in 1:3) {
    theta[k] ~ dnorm(0,0.00001)  # Vague normal on mean
    sigma[k] ~ dunif(0,100)      # Uniform on sigma
    tau0[k] <- 1/(sigma[k]*sigma[k])
  }
}

# Data:
list(N= 10, Y=c(1.2,2.4,1.3,1.3,0.0,1.0,1.8,0.8,4.6,1.4))


# Inits
list(theta = c(2,2,2), sigma = c(1,1,1), 
inv.lambda2=c(1,1,1,1,1,1,1,1,1,1),
lambda3=c(1,1,1,1,1,1,1,1,1,1)
)
list(theta = c(1,1,1), sigma = c(2,2,2), 
inv.lambda2=c(10,10,10,10,10,10,10,10,10,10),
lambda3=c(.1,.1,.1,.1,.1,.1,.1,.1,.1,.1)
)