branching

branching.R

pcr <- function(y,p,rounds) {
  for (i in 1:rounds) {
    y <- y + rbinom(1, y, p)
  }
  y
}

prob <- function(s,n,p) {
  if (s == 1) {
    if (n == 1) return(1)
    if (n != 1) return(0)
  }
  sum(sapply(0:ceiling(n/2), function(k) {
    prob(s-1,n-k,p) * choose(n-k, k) * p^k * (1-p)^(n-2*k)
  }))
}

p <- .7
rounds <- 5
hist(replicate(10000, pcr(1,p,rounds)), col="grey", 
     freq=FALSE, breaks=0:(2^rounds+1)-.5)
x <- 0:(2^rounds)
y <- sapply(x, function(n) prob(rounds+1,n,p))
points(x, y, col="blue", lwd=4, type="b")

# the mean of the distribution?
(1 + p)^rounds