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