######### MPT WinBUGS code to monitor the population mean

##################################################
# 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])
} # CHANGED BELOW TO NOT BE MEDIAN-MU ANYMORE AND UPDATE X[1,1]
# 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 1: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]))
g.prod[i]<-p[i]*m.g[i]
}
G.mean<-sum(g.prod[1:sets])*sigma+mu
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,
m.g=c(-0.0703,-0.0526,-0.0443,-0.0384,-0.0336,-0.0296,-0.0259,
-0.0226,-0.0195,-0.0166,-0.0139,-0.0112,-0.0086,-0.0061,
-0.0036,-0.0012,0.0012,0.0036,0.0061,0.0086,0.0112,0.0139,
0.0166,0.0195,0.0226,0.0259,0.0296,0.0336,0.0384,0.0443,0.0526,0.0703))