Poisson-binomial distribution - any existing Stan implementation?

Modeling
6

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)).