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

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:

• Unless you can also restrict the $p_i$ in a meaningful way, there's little one can say in general. This setup could be a discrete approximation to literally anything. – whuber Mar 28 '17 at 14:46
• Very good point. The data comes from real life measurements so it's unclear what exact distribution they follow. I've added charts showing the distributions of both the probabilities and weights. – Leon P Mar 30 '17 at 20:28
• Look into the saddlepoint approximation, search this site. I might be writing an answer along that line. – kjetil b halvorsen Mar 31 '17 at 0:10

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

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.

• The examples in the question itself clearly demonstrate this can be a bad approximation. There is no theory that generally suggests otherwise, because the Central Limit Theorem does not necessarily hold (it imposes specific conditions on the $w_i$ and $p_i$). – whuber Dec 9 '17 at 18:46
• @whuber, please explain how the examples given demonstrate that this can be a bad approximation. I don't get it. Thanks. – use_norm_approx Dec 9 '17 at 21:44
• I thought the answer given here (stats.stackexchange.com/questions/89254/…) implies everything is fine and dandy. – use_norm_approx Dec 9 '17 at 21:46
• Also I think this question is a duplicate of the other one. If not, I don't understand why. Thanks. – use_norm_approx Dec 9 '17 at 21:47
• The reference you found is a good one, and it says explicitly that you have to verify the conditions of the CLT (whether the Lindeberg or Lyapunov) before you can justify applying it. Due the great range of weights quoted in the present question, we have to believe this could be the kind of "weird edge case" mentioned by @Glen_b in his answer there. – whuber Dec 9 '17 at 22:53

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")
`