DataIn <- c(10,10,10,15,5)
Trans <- c(10,1,5,4)

mcmc <- function(datain,transin,pcut,ncap,burn,length){
# datain <- c(#no data, # early S, # early F, # final S, # final F)
# transin <- c(#Early S -> Final S, #Early S -> Final F, 
#                   Early F -> Final S, #Early F -> Final F) 
# pcut is threshold for Pr(pi > 0.5)
# ncap is maximum sample size
# burn is the mcmc burn in sample size
# length is the number of draws after burn
# output is pwin, which is (predict win @ N, predict win @ Cap)

  pwin <- numeric(2)

  pi <- 0.5
  gamma0 <- 0.5
  gamma1 <- 0.5

  for (i in 1:(burn+length)){
# Complete conditionals
    gamma0 <- rbeta(1,1+transin[3],1+transin[4])
    gamma1 <- rbeta(1,1+transin[1],1+transin[2])

    success <- datain[4]
    n <- sum(datain)
    success <- success + rbinom(1,datain[1],pi)
    success <- success + rbinom(1,datain[3],gamma0)
    success <- success + rbinom(1,datain[2],gamma1)

    pi <- rbeta(1,1+success,1+n-success)

# Predictive Probabilities: First win on current N
    probHalfN <- 1.0-pbeta(0.5,1+success,1+n-success)
    winTrialN <- ifelse(probHalfN > pcut,1,0)

# Predictive Probabilities: Second win at cap
    successCap <- success + rbinom(1,ncap-n,pi)
    probHalfCap <- 1.0-pbeta(0.5,1+successCap,1+ncap-successCap)
    winTrialCap <- ifelse(probHalfCap > pcut,1,0)

# Now record summaries of that past burn
    if (i > burn){
      pwin[1] <- pwin[1] + winTrialN
      pwin[2] <- pwin[2] + winTrialCap
    }
  }
  pwin <- pwin/length
  pwin
}

#> mcmc(DataIn,Trans,0.965,100,1000,10000)
#[1] 0.9292 0.9419