##############################################################
#
# The school coaching data - basic model
#
##############################################################

list(J=8,
  y=c(28, 8, -3, 7,  -1, 1, 18, 12),
  sigma.y = c(15, 10, 16, 11,  9, 11, 10, 18))

model {
  for (j in 1:J){
  y[j] ~ dnorm (theta[j], precis.y[j])   # note: N(mean,precis)
    theta[j] ~ dnorm (mu, precis.theta)  # where precis = 1/var
    precis.y[j] <- pow(sigma.y[j], -2)
  }
  mu ~ dnorm (0.0, 1.0E-6)
  precis.theta <- pow(tau, -2)
  tau ~ dunif (0, 1000)        # prior is p(m,t) = 1
}                              # rather than p(m,t^2) = IG(1,1)

##############################################################
#
# Darwin's cross-fertilization data: Distortion in the normal
# caused by outliers; t-modeling of outliers
#
# x[i] = difference in height between corn plants in pair i,
#        where one plant is cross-fertilized and the other is
#        self-fertilized.
#
# Our main interest is in mu = mean difference,and V = variance 
# of difference. 
#
##############################################################

# Model 1 (basic normal model):

model { 
  for (i in 1:15) { 
    x[i] ~ dnorm(mu,tau)  # note tau = precis
  }
  tau ~ dgamma(1,0.001)
  V <- 1/tau
  mu ~ dnorm(0,0.0000001)
}

# Model 1 Data: 

list(x=c(49,-67,8,16,6,23,28,41,14,29,56,24,75,60,-48))

# Model 1 Initial Values

list(mu=20,tau=1)

########################################################

# Model 2 (student t for x; estimate nu=df also!):

model { 
  for (i in 1:15) { 
    x[i] ~ dt(mu,tau,nu)
  }
  nu ~ dexp(k) I(2,)
  tau ~ dgamma(0.0001,0.0001)
  k ~ dunif(0.01,0.5)
  V <- 1/tau
  mu ~ dnorm(0,0.0000001)
}

# Model 2 Data: 

list(x=c(49,-67,8,16,6,23,28,41,14,29,56,24,75,60,-48))

# Model 2 Inits:

list(mu=0,tau=1,k=0.25)
list(mu=27,tau=0.001,k=0.11)
list(mu=46,tau=0.005,k=0.47)

########################################################

# Model 3 (same as model 2, but t = scale mixture of normals):

model { 
  for (i in 1:15) { 
    x[i] ~ dnorm(mu,tau[i])
    tau[i] <- lam[i]*Tau
    lam[i] ~ dgamma(nu,nu)
  }
  nu ~ dexp(k) I(2,)
  k ~ dunif(0.01,0.5)
  Tau ~ dgamma(0.0001,0.0001)
  V <- 1/Tau
  mu ~ dnorm(0,0.0000001)
}

# Model 3 Data: 

list(x=c(49,-67,8,16,6,23,28,41,14,29,56,24,75,60,-48))

# Model 3 Inits:

list(mu=0,Tau=1)
list(mu=24,Tau=0.001)
list(mu=47,Tau=0.002)

########################################################