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:

Values of P enter image description here enter image description here

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:

  • 4
    $\begingroup$ 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. $\endgroup$
    – whuber
    Commented Mar 28, 2017 at 14:46
  • $\begingroup$ 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. $\endgroup$
    – Leon P
    Commented Mar 30, 2017 at 20:28

4 Answers 4


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

  • $\begingroup$ It's a nice idea--but wouldn't you first have to establish that the distribution is unimodal (at least to a suitable approximation)? How would one go about doing that? Simulation? $\endgroup$
    – whuber
    Commented May 11, 2020 at 19:14
  • $\begingroup$ Simulations already shown in another answer is unimodal. But it would of course depend on distributions of parameters and weights. But this link seems to give some such approximations for bimodal distributions. (I will add some examples) $\endgroup$ Commented May 11, 2020 at 19:28
  • $\begingroup$ The "histogram of p" plot in the question is bimodal. Multimodality might be the norm: consider a situation with a few relatively large weights attached to variables with fairly different values of $p_i,$ to which are added many "noise"-like terms with small weights. $\endgroup$
    – whuber
    Commented May 11, 2020 at 20:13

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.

  • 3
    $\begingroup$ 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$). $\endgroup$
    – whuber
    Commented Dec 9, 2017 at 18:46
  • $\begingroup$ @whuber, please explain how the examples given demonstrate that this can be a bad approximation. I don't get it. Thanks. $\endgroup$ Commented Dec 9, 2017 at 21:44
  • $\begingroup$ I thought the answer given here (stats.stackexchange.com/questions/89254/…) implies everything is fine and dandy. $\endgroup$ Commented Dec 9, 2017 at 21:46
  • $\begingroup$ Also I think this question is a duplicate of the other one. If not, I don't understand why. Thanks. $\endgroup$ Commented Dec 9, 2017 at 21:47
  • 2
    $\begingroup$ 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 to 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. $\endgroup$
    – whuber
    Commented Dec 9, 2017 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
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)

enter image description here enter image description here

# define simulation function. `n_sim`` is number of simulations
sim_expr <- function(n_sim = 100000){
  # setup cluster
  cl <- makeCluster(7)
  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)
system.time(s1 <- sim_expr())
#R>   user  system elapsed 
#R>   0.11    0.03   14.22
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")

enter image description here


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.