model {
for (i in 1:20) {Y[i] ~ dgamma(4, 0.04)}
cum[1] <- Y[1]
for (i in 2:20) {
cum[i] <- cum[i - 1] + Y[i]
}
for (i in 1:20) {
cum.step[i] <- i*step(1000 - cum[i])
}
number <- ranked(cum.step[], 20) # maximum number in cum.step
check <- equals(cum.step[20], 0) # always 1 if I=20 big enough
}

node mean sd MC error 2.5% median 97.5% start sample
Y[1] 101.1 50.02 0.5408 28.51 93.23 222.2 1 10000
number 9.631 1.636 0.01503 7.0 10.0 13.0 1 10000