# How to efficiently estimate number of individuals with n+ successes from a series of bernoulli trials?

I have a situation where I need to estimate the number of persons exposed to a given event n or more times. For each person, I have an array of probabilities representing independent Bernoulli trials where a successful trial means 1 exposure of that individual to the event in question.

My data, for instance, looks something like this:

Person A: [0.34, 0.89, 0.01, 0.50]
Person B: [0.22, 0.45, 0.60]
Person C: [0.10, 0.10, 0.99, 0.32, 0.22]


To make it more concrete, if we pretend that I knew the outcome of the trials above as:

Person A: [0, 1, 1, 1]
Person B: [0, 1, 1]
Person C: [0, 0, 1, 0, 0]


...and the n I am calculating against is 2, then my answer is "2 people" because Person A and Person B were exposed to the event 2 or more times.

I know that the number of successful events for a series of independent Bernoulli trials (with different probabilities) is given by the Poisson Binomial Distribution, so for any given person, I could use standard methods for that distribution.

However, my understanding is that most methods of such have troubles or are untenable at a large number of probabilities. In my case, the array of probabilities for a given person may wildly vary (say 2 to 200 or more), and the number of persons can be prohibitively large (millions). Additionally I'm not just estimating "the number of events," but rather "the number of people with n+ events."

Is there a way to efficiently estimate (so for instance, not monte carlo) the above even if the estimate has an understood error propensity? For instance, at the moment I'm running some simulations where I simply sum up the probabilities per person and count up the persons where that is greater than my n which is sometimes close to the result a monte carlo would give, but I don't understand the error behavior yet or if there is a better method.

• It is presently unclear exactly what your data is ---i.e., wht you know for sure and what you are estimating.
– Ben
Commented Apr 24, 2022 at 0:30
• I'm afraid I don't follow? I have provided example data as well as an example estimation. If you can provide specifics on where the confusion arises, I'll endeavor to enhance the question.
– Josh
Commented Apr 24, 2022 at 1:00
• Happy to elaborate. You say that your data consists of probabilities in Bernoulli trials. It is unclear to me how you could observe this, so do you actually have known probabilities (and if so how?) or are you just observing some underlying binary data that has led you to those values as estimated probabilities?
– Ben
Commented Apr 24, 2022 at 2:05
• Arguably the probabilities I have are themselves indeed estimated ones, but my needs are content to not propagate whatever error comes from that for them moment. They arise from a separate, large, and rather convoluted process, but one well studied and that I trust. In either case, the above example data is indeed what I have (and need) to work with. I do not have the actual binary outcomes in the second/lower box which is just for demonstration of what I ultimately want to calculate.
– Josh
Commented Apr 24, 2022 at 17:10

## 1 Answer

I'm not convinced that the stipulated size for your data is too large for exact computation using the Poisson-binomial distribution. Nevertheless, since you are looking for some reasonable approximation method, you are probably going to have to fall back onto the normal approximation from the CLT (which ought to be slightly better then what you are doing now. Alternatively, you could consider using a hybrid method that I'll explain below.

Suppose we have a Poisson-binomial random variable $$X \sim \text{PoisBin}(\boldsymbol{\theta})$$ with probabilities $$\boldsymbol{\theta} = (\theta_1,...,\theta_k)$$ for $$k$$ independent trials. We can write this as $$X = \sum_i U_i$$ where $$U_1,...,U_k \sim \text{IID Bern}(\theta_i)$$, which allows us to easily compute the mean and variance:

\begin{align} \mathbb{E}(X) &= \mathbb{E} \bigg( \sum_{i=1}^k U_i \bigg) \\[6pt] &= \sum_{i=1}^k \mathbb{E} ( U_i ) \\[6pt] &= \sum_{i=1}^k \theta_i, \\[6pt] \mathbb{V}(X) &= \mathbb{V} \bigg( \sum_{i=1}^k U_i \bigg) \\[6pt] &= \sum_{i=1}^k \mathbb{V} ( U_i ) \\[6pt] &= \sum_{i=1}^k \theta_i(1-\theta_i). \\[6pt] \end{align}

Now, if $$k$$ is not too small, we can invoke the CLT to approximate the distribution by a normal distribution with the same mean and variance, giving the resulting approximation:

$$X \overset{\text{approx}}{\sim} \text{N}(\mu, \sigma^2) \quad \quad \quad \quad \quad \mu \equiv \sum_{i=1}^k \theta_i \quad \quad \quad \quad \quad \sigma^2 \equiv \sum_{i=1}^k \theta_i(1-\theta_i).$$

There are various ways you can obtain the specific approximating distribution you want to use (e.g., using continuous approximation with or without continuity correction, discrete approximation, etc.) and these will give you different formulae for the approximate probability of interest. I quite like to use the approximation obtained by taking the density at the exact discrete points and then normalising, which gives:

$$\mathbb{P}(X \geqslant r) \approx \frac{\sum_{x=r}^k \text{N}(x|\mu, \sigma^2)}{\sum_{x=0}^k \text{N}(x|\mu, \sigma^2)}.$$

(Note here that I use the notation $$r$$ instead of $$n$$ for the cut-off value of interest; it is bad practice to use the latter as notation, since it is the standard notation for the number of data points in the analysis.) It ought to be computationally feasible to compute this approximate probability for each of your rows. Each computation requires you to first compute $$\mu$$ and $$\sigma^2$$ from the above formulae and then compute the approximate probability

Implementation in R: You can implement this approximation in R quite simply using the existing facilities for the normal distribution. It is best to undertake computations in log-space for accuracy. Here is a function that will compute the approximate probability/log-probability for an input probability vector theta and input value n (with an option log to return the result as a log-probability instead of a probability).

#Create a function to compute the approximate probability
prob.approx <- function(theta, r, log = FALSE) {
k <- length(theta)
if (r > k) { if (log) { return(-Inf) } else { return(0) } }
MEAN <- sum(theta)
VAR  <- sum(theta*(1-theta))
LOGS <- dnorm(0:k, mean = MEAN, sd = sqrt(VAR), log = TRUE)
LOGP <- matrixStats::logSumExp(LOGS[(r+1):(k+1)]) - matrixStats::logSumExp(LOGS)
if (log) { LOGP } else { exp(LOGP) } }


We can experiment to see how fast this runs using the microbenchmark package. First I will set up some mock data where the number of elements in the probability vectors is Poisson distributed with a mean of 40 values. Then I will run a benchmarking test on a procedure that computes the approximating probability for $$n=1000$$ values (using a value for $$r$$ that is usually, but not always, below the number of trials).

#Set up some mock data (don't count this part in speed test)
set.seed(1)
n    <- 1000
LAM  <- 40
LENG <- rpois(n, lambda = LAM)
DATA <- vector(n, mode = 'list')
for (i in 1:n) { DATA[[i]] <- runif(LENG[i]) }

#Set up procedure
PROC <- function(data, r, log = FALSE) {
n   <- length(data)
OUT <- rep(0, n)
for (i in 1:n) { OUT[i] <-  prob.approx(data[[i]], r, log = log) }
OUT }

#Benchmark the procedure
library(microbenchmark)
set.seed(1)
microbenchmark(PROC(DATA, 30))

Unit: milliseconds
expr     min      lq     mean   median      uq     max neval
PROC(DATA, 30) 12.0072 12.2244 12.95248 12.27545 12.6327 26.7406   100


As can be seen from the output of this benchmarking test, it took about 13 milliseconds to compute the approximating probabilities for $$n=1000$$ data values. Consequently, it should be possible to compute for $$n=10^6$$ values in about 13 seconds. While I have not used quite the same spread of vector lengths as what you describe in your question, this method ought to be computationally feasible in your problem.

• I'm running a couple of sims to see if I can get this to play out, but this looks good so far and has better results than my naive approx (I think...). One question though, in your R code you slice the pdf results at: [(r+1):(k+1)] whereas I would expect the slice to be [r:k] or maybe [r+1:k] given the formula you posted. Can you elaborate? I'm not sure if this is me misunderstanding the statistics or a difference of flavor of Python background (in Python the first element is at index = 0)?
– Josh
Commented Apr 25, 2022 at 15:09
• Also yes, I think you're right one run of this may be sensible for an exact approximation via the typical ways of calculating a Poisson Binomial. I forgot to add in my question though that I'm likely to have to do this calc many many times (with possibly millions of persons and varying probability array lengths) to explore some different possibilities, but otherwise yep would be better just to wait for a solo calc of it!
– Josh
Commented Apr 25, 2022 at 15:13
• The vector LOGS has k+1 elements starting at a value of zero, so the elements (r+1):(k+1) of that vector correspond to the values $r,...,k$.
– Ben
Commented Apr 25, 2022 at 22:18
• Odd, In Python I'm just slicing as [r:k] as I don't seem to have that extra element (k+1 doesn't exist). Granted I'm skipping the log bit and using scipy's norm.pdf(), but it seems to work out well. Sims also verify this is way better than my current methods; however, as expected it struggles with low k, albeit I might be able to use a hybrid approach (calc directly for low k and norm approx for high k). Either way this is great and gets me enough to move forward!
– Josh
Commented Apr 26, 2022 at 15:05