########## the code that doesn't monitor \int y G ( dy );
# instead it monitors the population median

##################################################
# 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 )
left <- mean( y[ ] ) - 3.25 * sd( y[ ] )
right <- mean( y[ ] ) + 3.25 * sd( y[ ] )
}
##################################################
# starting values
##################################################
list( mu = 404.59, sigma = 6.458053, c = 2 )
##################################################
# data; n=100 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(387.9325, 390.5564, 391.9152, 392.8727, 393.6261, 394.2547,
  394.7985, 395.2808, 395.7162, 396.1147, 396.4834, 396.8275, 397.1509,
  397.4567, 397.7473, 398.0247, 398.2906, 398.5462, 398.7926, 399.0310,
  399.2620, 399.4864, 399.7048, 399.9178, 400.1259, 400.3294, 400.5288,
  400.7244, 400.9165, 401.1054, 401.2914, 401.4747, 401.6556, 401.8342, 
  402.0107, 402.1853, 402.3581, 402.5294, 402.6993, 402.8678, 403.0352,
  403.2016, 403.3670, 403.5316, 403.6956, 403.8590, 404.0219, 404.1845,
  404.3468, 404.5089, 404.6711, 404.8332, 404.9955, 405.1581, 405.3210, 
  405.4844, 405.6484, 405.8130, 405.9784, 406.1448, 406.3122, 406.4807, 
  406.6506, 406.8219, 406.9947, 407.1693, 407.3458, 407.5244, 407.7053, 
  407.8886, 408.0746, 408.2635, 408.4556, 408.6512, 408.8506, 409.0541, 
  409.2622, 409.4752, 409.6936, 409.9180, 410.1490, 410.3874, 410.6338, 
  410.8894, 411.1553, 411.4327, 411.7233, 412.0291, 412.3525, 412.6966, 
  413.0653, 413.4638, 413.8992, 414.3815, 414.9253, 415.5539, 416.3073, 
  417.2648, 418.6236, 421.2475),
n=100, Jt=5,sets=32,
j.int=c(1,2,4,8,16))