Probability that sum of binary variables is even

Question

Let $S_i \in \{0,1\}$, $i=1,\dots,N$ be $N$ independent random binary variables, each taking the value 1 with probability $0 \le p_i \le 1$ (and the value 0 with probability $1-p_i$).

I am interested in the sum,

$$X = \sum_{i=1}^N S_i$$

Because the $p_i$ can be different, the distribution of $X$ can be complicated (https://en.wikipedia.org/wiki/Poisson_binomial_distribution).

I'm not interested in the full distribution of $X$, but only on the probability that $X$ is even (divisible by 2). Can this probability be computed in a simple (or efficient) manner?

Answer 1

You should be able to compute this quickly through a very simple dynamic programming approach.

Let $q_i$ denote the probability that $S_1+\dots+S_i$ is even.

Then $q_1=1-p_1$.

Going from $i-1$ to $i$, you either have $S_1+\dots+S_{i-1}$ odd with probability $1-q_{i-1}$, then $S_1+\dots+S_i$ will be even if $S_i=1$, for a total probability of $(1-q_{i-1})p_i$, or the other way around. Overall,

$$ q_i = (1-q_{i-1})p_i + q_{i-1}(1-p_i). $$

Just calculate $q_N$ by iterating over $i$, and there you are.

Answer 2

There's a nice approach that uses generating functions. If $G(z) = \mathbf{E} z^X$ is the probability generating function of a non-negative random variable $X$, then the probability that $X$ is even is $(1 + G(-1))/2$, as per this discussion: https://math.stackexchange.com/questions/1634021/even-value-in-a-r-v-using-generating-functions.

But we have that \begin{align*} G(z) = \mathbf{E} z^{\sum_{i=1}^N S_i} = \prod_{i=1}^N \mathbf{E} z^{S_i} = \prod_{i=1}^N (p_i z + (1 - p_i)), \end{align*} from which we can calculate \begin{align*} \mathbf{P}(X \text{ is even}) = \frac{1}{2}(1 + G(-1)) = \frac{1}{2}\Bigl(1 + \prod_{i=1}^N (1 - 2p_i)\Bigr). \end{align*}

A fun observation from the above is that if even one of the $p_i$ are equal to $1/2$, then the probability of $X$ being even is $1/2$, regardless of the other probabilities. This is because $S_i$ changing from $0$ to $1$ (or $1$ to $0$) switches the parity of the whole sum, and that that switch happens with probability $1/2$ independent of the rest of the sum.

Answer 3

There is an elementary solution requiring almost no work to derive.

---

(I am going to switch the notation because "$S$" is such a handy mnemonic for a sum.)

Let $f_k$ be the chance that the sum $S_k = X_1+X_2+\cdots+X_k$ is even. The event "$S_{k+1}$ is even" can be expressed as the disjoint union of "$S_k$ is odd and $X_{k+1}$ is odd" and "$S_k$ is even and $X_{k+1}$ is even." Thus, their probabilities add; and since $X_{k+1}$ is assumed independent of $S_k,$ the chances of these two component events are computed via multiplication to obtain

$$f_{k+1} = (1 - f_k)p_{k+1} + f_k(1 - p_{k+1}) = p_{k+1} + (1 - 2p_{k+1})f_k.$$

This is the relation reported by Stephan Kolassa elsewhere in this thread. It is a first-order recursion beginning with $f_0 = 1$ (an empty sum is defined to be zero, which is certainly even).

There are many ways to solve it. One of the simplest is to write $q_k = 1 - 2p_k$ and express the recursion in terms of $a_k = 2f_k - 1$ as

$$a_{k+1} = q_{k+1}a_k;\quad a_0 = 2(1) - 1 = 1.$$

This series obviously has the unique solution $a_k = (1)(q_1)(q_2)\cdots(q_k),$ whence (returning to the original notation in terms of $f_k$ and $p_k$)

$$f_k = (1 + a_k)/2 = \frac{1}{2}\left(1 + \prod_{i=1}^k (1 - 2p_i)\right),$$

exactly the formula beautifully obtained by Damian Pavlyshyn elsewhere in this thread.

Answer 4

Here is a solution to the generalised problem using other modulo

The other answers here are wonderful, but I can't help but scratch the itch to generalise your problem. Consider the generalised problem where we have a sequence $X_1,X_2,X_3,...$ with independent binary values having $X_i \sim \text{Bern}(p_i)$. Now define the modular-remainders of the sums of interest as:

$$M_{n,k} \equiv\bigg( \sum_{i=1}^n X_i \bigg) \text{ mod } k$$

for $n=1,2,3,...$ and $k=2,3,4,...$. We want to compute the probability distribution for this statistic, with probabilities denoted by:

$$\psi_{n,k}(m) \equiv \mathbb{P}(M_{n,k} = m) \quad \quad \quad \quad \quad \text{for } m =0,...,k-1.$$

Your question is seeking the probability $\psi_{n,2}(0)$, which occurs in the case $k=2$ where the value $M_{n,2}$ is an indicator that the sum of the first $n$ values in the sequence is even.

---

Recursive computation of the probability of interest: In this general case we can use an extension of the recursive algorithm given in the wonderful answer by Stephen Kolassa. To do this, we can use the base equation:

$$\psi_{0,k}(m) = \mathbb{I}(m=0).$$

and the recursive equation:

$$\psi_{n+1,k}(m) = (1-p_{n+1}) \cdot \psi_{n,k}(m) + p_{n+1} \cdot \psi_{n,k} ((m-1) \text{ mod } k).$$

You can find a similar generalisation relating to an associated problem in this related answer. It is possible to obtain a closed-form solution for the probability function of interest (as opposed to a recursive formula) using the method described in the linked question. Here we will proceed using the recursive computation, which allows efficient computation of matrices of probabilities in cases where we are interested in all sample sizes up to some upper bound.

---

Implementation: We can implement this recursive method in R using a function that takes n, k and p as inputs and generating the vector of probabilities over all relevant values of $m$. (Note: In the function below, for the input p we can take this as either a vector of probabilities that is at least $n$ elements long, or a function that generates a probability for any positive integer input to allow a full sequence of probabilities.) The function defaults to the binary case that is of interest to you but it allows different values of k to be set by the user.

modular.prob <- function(n, p, k=2, keep.all = FALSE) {
  
  #Check inputs n and k
  if (!is.numeric(n))      stop('Error: Input n should be a number')
  if (!is.numeric(k))      stop('Error: Input k should be a number')
  if (length(n) != 1)      stop('Error: Input n should be a single number')
  if (length(k) != 1)      stop('Error: Input k should be a single number')
  if (n != as.integer(n))  stop('Error: Input n should be an integer')
  if (k != as.integer(k))  stop('Error: Input k should be an integer')
  if (min(n < 0))          stop('Error: Input n should be a non-negative integer')
  if (min(k < 1))          stop('Error: Input n should be a positive integer')
  
  #Check input p
  if (is.function(p)) {
    if (length(formals(p)) != 1) {
      stop('Error: If input p is a function, it should only have a single argument') }
    pp <- rep(0, n)
    for (i in 1:n) { pp[i] <- p(i) }
  } else {
    if (!is.numeric(p))    stop('Error: Input p should be a numeric vector')
    if (length(p) < n)     stop('Error: Input p should have at least n elements')
    pp <- p[1:n]
    if (min(pp) < 0)       stop('Error: Input p should have values between zero and one')
    if (max(pp) > 1)       stop('Error: Input p should have values between zero and one') }
  
  #Check input keep.all
  if (!is.logical(keep.all)) stop('Error: Input keep.all should be a logical value')
  if (length(keep.all) != 1) stop('Error: Input keep.all should be a single logical value')
  
  #Set probability matrix
  PROBS <- matrix(0, nrow = n+1, ncol = k)
  rownames(PROBS) <- sprintf('n[%s]', 0:n)
  colnames(PROBS) <- sprintf('m[%s]', 1:k-1)
  if (k == 2) { colnames(PROBS) <- c('Even', 'Odd') }
  
  #Compute probabilities
  #Trivial case where k = 1 is treated without recursion
  PROBS[1, 1] <- 1
  if (k == 1) { PROBS[, 1] <- 1 }
  if (k > 1) { 
  for (i in 1:n) {
    PP <- PROBS[i, ]
    PROBS[i+1, ] <- (1-pp[i])*PP + pp[i]*PP[c(k, 1:(k-1))] } }
   
   #Give output
   if (keep.all) { PROBS } else { PROBS[n+1, ] } }

Here is an example implementing the function with a specified probability vector with $n=20$ elements. We first compute the probabilities for even/odd sums using $k=2$ and we then compute the probabilities for modular remainders with $k=3$. In each case we get a matrix of probabilities for each trial and each modular remainder.

#Set probability vector (with twenty elements)
p <- c(0.1, 0.2, 0.3, 0.2, 0.1, 0.6, 0.2, 0.2, 0.5, 0.1, 
       0.9, 0.2, 0.2, 0.1, 0.4, 0.5, 0.6, 0.4, 0.2, 0.8)

#Compute probabilities of even/odd sums
modular.prob(n = 20, p = p, keep.all = TRUE)

           Even       Odd
n[0]  1.0000000 0.0000000
n[1]  0.9000000 0.1000000
n[2]  0.7400000 0.2600000
n[3]  0.5960000 0.4040000
n[4]  0.5576000 0.4424000
n[5]  0.5460800 0.4539200
n[6]  0.4907840 0.5092160
n[7]  0.4944704 0.5055296
n[8]  0.4966822 0.5033178
n[9]  0.5000000 0.5000000
n[10] 0.5000000 0.5000000
n[11] 0.5000000 0.5000000
n[12] 0.5000000 0.5000000
n[13] 0.5000000 0.5000000
n[14] 0.5000000 0.5000000
n[15] 0.5000000 0.5000000
n[16] 0.5000000 0.5000000
n[17] 0.5000000 0.5000000
n[18] 0.5000000 0.5000000
n[19] 0.5000000 0.5000000
n[20] 0.5000000 0.5000000

#Compute probabilities for ternary modular remainders
modular.prob(n = 20, p = p, k = 3, keep.all = TRUE)

           m[0]      m[1]      m[2]
n[0]  1.0000000 0.0000000 0.0000000
n[1]  0.9000000 0.1000000 0.0000000
n[2]  0.7200000 0.2600000 0.0200000
n[3]  0.5100000 0.3980000 0.0920000
n[4]  0.4264000 0.4204000 0.1532000
n[5]  0.3990800 0.4210000 0.1799200
n[6]  0.2675840 0.4078480 0.3245680
n[7]  0.2789808 0.3797952 0.3412240
n[8]  0.2914294 0.3596323 0.3489382
n[9]  0.3201838 0.3255309 0.3542853
n[10] 0.3235940 0.3249962 0.3514098
n[11] 0.3486283 0.3237342 0.3276375
n[12] 0.3444301 0.3287130 0.3268569
n[13] 0.3409155 0.3318564 0.3272281
n[14] 0.3395467 0.3327623 0.3276909
n[15] 0.3348044 0.3354761 0.3297195
n[16] 0.3322620 0.3351403 0.3325978
n[17] 0.3324635 0.3334133 0.3341233
n[18] 0.3331274 0.3330333 0.3338393
n[19] 0.3332698 0.3330522 0.3336781
n[20] 0.3335964 0.3332262 0.3331773

As exhibited in the outputs above, typically we will see that $M_{n,k}$ converges in distribution to the uniform distribution as $n \rightarrow \infty$. This can be proved formally using the CLT - it occurs due to the fact that the standardised sum of the binary values converges in distribution to a continuous distribution (the normal distribution) and so the modular-remainder over shorter and shorter intervals converges to uniform probabilities.