####################
#  WinBUGS code for DIC comparison of normal scale mixture models
#  See Carlin and Louis (2008), Example 4.8
#################### 


## Compare normal scale mix and "direct fit" of t_2 and DE error structures:
model
{
	for(i in 1:N) {
		# Duplicate the data
		Y1[i] <- Y[i]
		Y2[i] <- Y[i]
		Y3[i] <- Y[i]
		Y4[i] <- Y[i]
		Y5[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
		lambda3[i] ~ dexp(0.5);           # M3 = e_i ~ DE(0,tau3)
		
		# Direct fit the t_2 and DE
		Y4[i] ~ dt(theta[4],tau0[4],2);
		Y5[i] ~ ddexp(theta[5],tau0[5]);
	}
	# Priors
	for (k in 1:5) {                   
		theta[k] ~ dnorm(0,0.00001)      # Vague normal prior on grand slope
		sigma[k] ~ dunif(0,100)          # Uniform prior on sigma
		tau0[k] <- 1/(sigma[k]*sigma[k])
	}
}


## Sleep 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,2,2), sigma = c(1,1,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,1,1), sigma = c(2,2,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)
)