model {
for(i in 1:N) {
y[i] ~ dbern(theta)
y.rep[i] ~ dbern(theta)
}
theta ~ dunif(0,1)
switch.obs[1] <- 0
switch.rep[1] <- 0
for (i in 2:N) {
switch.obs[i] <- 1 - equals(y[i-1], y[i])
switch.rep[i] <- 1 - equals(y.rep[i-1], y.rep[i])
}
s.obs <- sum(switch.obs[])
s.rep <- sum(switch.rep[])
P.switch <- step(s.obs-s.rep-0.5)
+ 0.5*equals(s.obs, s.rep)
run.obs[1] <- 1
run.rep[1] <- 1
max.run.obs[1] <- 1
max.run.rep[1] <- 1
for (i in 2:N) {
run.obs[i] <- 1 + run.obs[i-1]
* equals(y[i-1], y[i])
run.rep[i] <- 1 + run.rep[i-1]
* equals(y.rep[i-1], y.rep[i])
max.run.obs[i] <- max(max.run.obs[i-1], run.obs[i])
max.run.rep[i] <- max(max.run.rep[i-1], run.rep[i])
}
P.run <- step(max.run.obs[N]-max.run.rep[N]-0.5)
+ 0.5*equals(max.run.obs[N], max.run.rep[N])
}
Data:
list(
N = 30,
y=c(0,0,0,1,0,1,0,0,0,1,0,0,1,0,1,0,0,0,0,1,1,0,0,0,1,0,0,1,0,1))
node mean sd MC error 2.5% median 97.5% start sample
P.run 0.1153 0.256 8.294E-4 0.0 0.0 1.0 1 100000
P.switch 0.9094 0.2642 8.941E-4 0.0 1.0 1.0 1 100000
max.run.rep[30] 6.937 3.047 0.008692 3.0 6.0 15.0 1 100000
s.rep 12.71 3.259 0.01033 6.0 13.0 19.0 1 100000
theta 0.3442 0.08273 2.597E-4 0.192 0.3408 0.5151 1 100000