R/poisson_binomial.R

In vlulla/vlutils: Common utilities that I find useful

poisson_binomial <-
function(theta) {
  N <- length(theta)
  if (N == 0) return(c(1))
  alpha <- matrix(-Inf, N+1, N+1)
  alpha[1,1] <- 1
  for(n in 1:N) {
    tot = 0
    alpha[n+1, tot+1] = alpha[n,tot+1]*(1-theta[n])
    if (n > 1) {
      for (tot in 1:(n-1)) {
        alpha[n+1,tot+1] <- alpha[n,tot] * theta[n] + alpha[n, tot+1]*(1-theta[n])
      }
    }
    tot = n
    alpha[n+1,tot+1] = alpha[n,tot]*theta[n]
  }
  alpha[N+1,1:(N+1)]
}