##################################################
# MPT code, uses "ones" trick
# c ~ gamma(5,1), mu ~ dnorm(0,0.001)
# sigma ~ dgamma(0.01,0.01) 
##################################################
model{
 for(i in 1:n){
  ones[i] <- 1
  z[i] <- (y[i]-mu)/sigma
  ind[i] <- trunc(sets*phi(z[i])+1)
  like[i] <- coi*exp(-0.5*z[i]*z[i])*p[ind[i]]/sigma
  # E{cpo.inv[i]|data}=1/CPO[i]
  cpo.inv[i] <- 1/(like[i]*sets)
  ones[i] ~ dbern(like[i])
 }
 x[1,1] <- 0.5; x[1,2] <- 0.5
   for(k in 1:j.int[Jt]){
    m[1,k] <- log(x[1,1])
    m[1,j.int[Jt]+k] <- log(x[1,2])
  }
 for(i in 2:Jt){
  par.lev[i] <- c*i*i
  for(j in 1:j.int[i]){
   x[i,2*j-1] ~ dbeta(par.lev[i],par.lev[i]); x[i,2*j] <- 1-x[i,2*j-1]
   for(k in 1:j.int[Jt-i+1]){
    m[i,(2*j-2)*j.int[Jt-i+1]+k] <- log(x[i,2*j-1])
    m[i,(2*j-1)*j.int[Jt-i+1]+k] <- log(x[i,2*j])
  }
 }
}
for(i in 1:sets){
 p[i] <- exp(sum(m[1:Jt,i]))
 }
for(i in 1:200){
 grid[i]<-left+(right-left)*i/201
 z.f[i] <- (grid[i]-mu)/sigma
 ind.f[i] <- trunc(sets*phi(z.f[i])+1)
 # f[i] is predictive density over grid[i]
 f[i] <- sets*coi*exp(-0.5*z.f[i]*z.f[i])*p[ind.f[i]]/sigma
 }
mu ~ dnorm(0,0.001)
sigma ~ dgamma(0.01,0.01) 
coi <- 0.3989422804
c ~ dgamma(5,1)
}

##################################################
# starting values
##################################################
list(mu=21,sigma=6,c=2)

##################################################
# data
# n=82 is sample size
# Jt=5 is level of finite Polya tree
# j.int is number of independent beta distributions at each level
# 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, Jt=5,sets=32,
j.int=c(1,2,4,8,16),left=0,right=40)