##################################################
# Dirichlet-multinomial allocation model for density estimation
# Approximates a DP(alpha G0)
# Fixed alpha (makes updating much faster)
# Common precision tau for each component
##################################################
model{
 for(i in 1:n){
  y[i] ~ dnorm(mu[z[i]],tau)
  z[i] ~ dcat(prob[])
  for(j in 1:K){
   f.d[i,j] <- prob[j]*coi*exp(-0.5*pow(y[i]-mu[z[i]],2)*tau)*sqrt(tau)
  }
  # E{cpo.inv[i]|data}=1/CPO[i]
  cpo.inv[i] <- 1/sum(f.d[i,1:K])
 }
 prob[1:K] ~ ddirch(prob.par[])
 for(i in 1:K){
  repeat[i] <- step(1-dist[i])
  prob.par[i] <- alpha/K
  mu[i] ~ dnorm(mu.mean,mu.prec)
 }
 # total is number of active "used" components
 total <- sum(repeat[])
 dist[1] <- 1
 for(i in 2:K){
  dist[i] <- sum(count[i,1:i])
  count[i,i] <- 1
  for(j in 1:i-1){
   count[i,j] <- equals(z[i],z[j])
  }
 }
 tau ~ dgamma(0.001,0.001)
 mu.mean ~ dnorm(0,0.001)
 mu.prec ~ dgamma(0.001,0.001)
 for(i in 1:200){
  grid[i]<-left+(right-left)*i/201
  for(j in 1:K){
   f.p[i,j] <- prob[j]*coi*exp(-0.5*pow(grid[i]-mu[j],2)*tau)*sqrt(tau)
  }
  # f[i] is predictive density over grid[i]
  f[i] <- sum(f.p[i,1:K])
  }
coi <- 0.3989422804
}

##################################################
# starting values
##################################################
list(mu.mean=21,mu.prec=0.03,tau=1.7)

##################################################
# data
# n=82 is sample size
# K is number of components in approximation
# alpha is Dirichlet process precision
# left and right are interval endpoints to evaluate predictive density
##################################################
list(y=c(9.172,9.35,9.483,9.558,9.775,10.227,10.406,16.084,
16.17,18.419,18.552,18.6,18.927,19.052,19.07,19.33,19.343,
19.349,19.44,19.473,19.529,19.541,19.547,19.663,19.846,19.856,
19.863,19.914,19.918,19.973,19.989,20.166,20.175,20.179,20.196,20.215,20.221,20.415,20.629,20.795,20.821,20.846,20.875,20.986,21.137,21.492,21.701,21.814,21.921,21.96,22.185,22.209,22.242,
22.249,22.314,22.374,22.495,22.746,22.747,22.888,22.914,23.206,23.241,23.263,23.484,23.538,23.542,23.666,23.706,23.711,24.129,24.285,24.289,24.366,24.717,24.99,25.633,26.96,26.995,32.065,32.789,34.279), 
n=82,K=50,alpha=1,left=0,right=40)