model{
   #Likelihood:
   for(i in 1:N){
     Y[i]~dnorm(mu[i],tau)
     mu[i]<-int+inprod(x[i,],beta[])
   }

   #Prior:
   tau~dgamma(0.1,0.1)
   int~dnorm(0,0.001)
   for(j in 1:p){
      beta[j]~dnorm(0,0.1)
   }

   #Diagnostics
   for(i in 1:N){
     Yrep[i]~dnorm(mu[i],tau)
     E[i]<-Y[i]-mu[i]
     Erep[i]<-Yrep[i]-mu[i]
     DY[i]<-step(Yrep[i]-Y[i])
   }
   Min_Y_samp<-ranked(Y[],1)
   Max_Y_samp<-ranked(Y[],N)
   Min_Y_rep <-ranked(Yrep[],1)
   Max_Y_rep <-ranked(Yrep[],N)
   D[1]<-step(Min_Y_rep-Min_Y_samp)
   D[2]<-step(Max_Y_rep-Max_Y_samp)
   D[3]<-step((Max_Y_rep-Min_Y_rep)-(Max_Y_samp-Min_Y_samp))

   Min_E_samp<-ranked(E[],1)
   Max_E_samp<-ranked(E[],N)
   Min_E_rep <-ranked(Erep[],1)
   Max_E_rep <-ranked(Erep[],N)
   D[4]<-step(Min_E_rep-Min_E_samp)
   D[5]<-step(Max_E_rep-Max_E_samp)
   D[6]<-step((Max_E_rep-Min_E_rep)-(Max_E_samp-Min_E_samp))

}

#Initial values:
list(beta=c(0,0,0),tau=1,int=15)

#Stacks data:
list(p = 3,   N = 21,  
	Y = c(42, 37, 37, 28, 18, 18, 19, 20, 15, 14, 14, 13, 11, 12, 8, 7, 8, 8, 9, 15, 15),
	x = structure(.Data = c(80, 27, 89,
				80, 27, 88,
				75, 25, 90,
				62, 24, 87,
				62, 22, 87,
				62, 23, 87,
				62, 24, 93,
				62, 24, 93,
				58, 23, 87,
				58, 18, 80,
				58, 18, 89,
				58, 17, 88,
				58, 18, 82,
				58, 19, 93,
				50, 18, 89,
				50, 18, 86,
				50, 19, 72,
				50, 19, 79,
				50, 20, 80,
				56, 20, 82,
				70, 20, 91),
 .Dim = c(21, 3)))