What is the CDF of the sum of weighted Bernoulli random variables?

Question

Let's say we have a random variable $Y$ defined as the sum of $N$ Bernoulli variables $X_i$, each with a different, success probability $p_i$ and a different (fixed) weight $w_i$. The weights are positive and between ~0.001-1,000 with a handful of outliers. If it would make the problem more tractable, it would be OK to discard both the very large and very small outliers.

Formally, $Y = \sum X_i w_i$

Where $\Pr(X_i=1)=p_i$ and $\Pr(X_i=0)=1-p_i$.

$N$ can range from ~10-10,000. I would like to quickly compute a approximation to $\Pr(Y<=k)$ (where $k$ is given).

The data I'm handling comes from real life bets, each with implied odds $p_i$ and wager $w_i$. Here are plots showing the rough distribution of the probabilities and weights:

One way to compute the approximation is through Monte Carlo simulation, but that can get slow for large values of $N$. One can also use the Central Limit Theorem to compute this for large $N$, but the accuracy is quite poor for small $N$ (or if there are a handful of large weights). Finally, there are ways to improve on the CLT approach by using asymptoptic expansions, e.g.:

Volkova, A. Y. (1996). A refinement of the central limit theorem for sums of independent random indicators. Theory of Probability and its Applications 40, 791-794.

However, as far as I've seen, the refined approximations only specify how to compute $Y$ if all weights $W_i$ are equal to 1 (a standard Poisson Binomial distribution). Is there a way to compute an answer that works for small and large $N$ with a smaller error than a pure CLT approach without having to resort to Monte Carlo Simulation? In other words, is it possible to extend the improved closed-form approximations to allow weights not equal to 1?

Some related questions that don't quite answer mine:

• CLT can be used for weighted sum of different Bernoulli variables? (Pure CLT is too inaccurate)
• How to find the distribution of the weighted sum of independent Bernoulli random variables for positive non-integer weights (No answer)
• https://math.stackexchange.com/questions/1546366/distribution-of-weighted-sum-of-bernoulli-rvs (The bounds mentioned in [Raghavan, 1988] are very, very wide for small N or when computing the CDF near the mean)

Answer 1

One possibility is to use the saddlepoint approximation. For that we need the mgf (moment generating function) and its logarithm the cgf (cumulant generating function.) The mgf of a Bernoulli variable with parameter $p_i$ is $$ \DeclareMathOperator{\E}{\mathbb{E}} M_i(t)=\E e^{t X_i}=(1-p_i)+p_i e^t $$ Then $$ \E e^{tY}=\prod_i \E e^{t w_i}=\prod_i M_i(t w_i) $$ and the cgf is $$ K(t)=\sum_{i=1}^n \log\left( 1-p_i+e^{t w_i} \right) $$ Then we can develop the saddlepoint approximation following How does saddlepoint approximation work?.

Answer 2

Use the normal approximation with

$$ \mu = \sum_{i=1}^n p_i w_i $$

and variance

$$ \sigma^2 = \sum_{i=1}^n w_{i}^2 p_i (1-p_i) $$

This should be pretty accurate, especially if $n$ is large.

Answer 3

Though you do ask for solution without Monte Carlo simulation, then I will highlight that it is not too computationally expensive in this case particularly if you use parallel computing. Here is an R example

#####
# simulate data
set.seed(91874947)
N  <- 10000 
ps <- ifelse(.4 > runif(N), rbeta(N, 5, 8), rbeta(N, 14, 5)) 
w  <- exp(rnorm(N, -1, .66) + (.2 > runif(N)) * rnorm(N, 1.4, .2))

hist(    ps, breaks = 50)
hist(log(w), breaks = 50)

#####
# define simulation function. `n_sim`` is number of simulations
sim_expr <- function(n_sim = 100000){
  # setup cluster
  require(parallel)
  cl <- makeCluster(7)
  on.exit(stopCluster(cl))
  clusterSetRNGStream(cl)
  clusterExport(cl, c("ps", "w", "N"))

  # run simulation
  parSapply(cl, 1:n_sim, function(...){
    y <- ps > runif(N)
    sum(w * y) 
  })
}

# run simulations and check run time. See the `elapsed` time (time in seconds)
set.seed(37219838)
system.time(s1 <- sim_expr())
#R>   user  system elapsed 
#R>   0.11    0.03   14.22
set.seed(4382482)
system.time(s2 <- sim_expr())
#R>   user  system elapsed 
#R>   0.08    0.00   14.11

# yields almost the same density estimate
plot (density(s1))
lines(density(s2), col = "red")

Answer 4

My intuition would suggest the following variation of Monte Carlo (although I cannot think of a proof at the minute). Start from a given realization of $X$, compute its $Y$ and check if less than $k$. Then choose an index $i$ with probability $w_i/(\Sigma_j w_j)$ (easy and efficient to do by just sampling a uniform variate and use bisection to find its place in the vector cumsum of $w_i/(\Sigma_j w_j)$ that you compute once at the beginning), resample $X_i$ from its distribution and if it has changed update $Y$ (less than one addition on average). Check again if less than $k$. Repeat this approach for several iterations and several starting points (akin to MCMC chains) keeping the tally of the times $Y$ was less than $k$. The approach is similar to Monte Carlo but selectively updates more often elements $X_i$ that tend to have a larger effect on $Y$.