Here’s a recursive version. The boundary conditions are somehow trickier than I was expecting:
dpois_binom <- function(y, theta) {
N <- length(theta);
if (N == 0) return(ifelse(y == 0, 1, 0));
if (N == 1) return(ifelse(y == 1, theta[1], 1 - theta[1]));
if (y == 0) return(prod(1 - theta));
# y >= 1, N >= 2
tail <- theta[2:N];
return((1 - theta[1]) * dpois_binom(y, tail) + theta[1] * dpois_binom(y - 1, tail));
}
This one just calculates a single value. But it’s still exponential in max(y, size(theta))
.