Probabilidad de que la suma de variables binarias sea par

Question

Sea $S_i \in \{0,1\}$, $i=1,\dots,N$ variables binarias aleatorias e independientes, cada una tomando el valor 1 con probabilidad $0 \le p_i \le 1$ (y el valor 0 con probabilidad $1-p_i$).

Me interesa la suma,

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

Debido a que los $p_i$ pueden ser diferentes, la distribución de $X$ puede ser complicada (https://en.wikipedia.org/wiki/Poisson_binomial_distribution).

No me interesa la distribución completa de $X$, sino solo la probabilidad de que $X$ sea par (divisible por 2). ¿Se puede calcular esta probabilidad de manera simple (o eficiente)?

Answer 1

Hay un enfoque interesante que utiliza funciones generadoras. Si $G(z) = \mathbf{E} z^X$ es la función generadora de probabilidad de una variable aleatoria no negativa $X$, entonces la probabilidad de que $X$ sea par es $(1 + G(-1))/2$, según esta discusión: https://math.stackexchange.com/questions/1634021/even-value-in-a-r-v-using-generating-functions.

Pero tenemos que \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*} desde donde podemos calcular \begin{align*} \mathbf{P}(X \text{ es par}) = \frac{1}{2}(1 + G(-1)) = \frac{1}{2}\Bigl(1 + \prod_{i=1}^N (1 - 2p_i)\Bigr). \end{align*}

Una observación interesante de lo anterior es que si alguna de las $p_i$ es igual a $1/2$, entonces la probabilidad de que $X$ sea par es $1/2$, independientemente de las otras probabilidades. Esto se debe a que el cambio de $S_i$ de $0$ a $1$ (o de $1$ a $0$) cambia la paridad de toda la suma, y ese cambio ocurre con probabilidad $1/2$ independientemente del resto de la suma.

Answer 2

Hay una solución elemental que requiere casi ningún trabajo para derivar.

---

(Voy a cambiar la notación porque "$S$" es un mnemotécnico tan útil para una suma.)

Sea $f_k$ la probabilidad de que la suma $S_k = X_1+X_2+\cdots+X_k$ sea par. El evento "$S_{k+1}$ es par" se puede expresar como la unión disjunta de "$S_k$ es impar y $X_{k+1}$ es impar" y "$S_k$ es par y $X_{k+1}$ es par." Por lo tanto, sus probabilidades se suman; y dado que se asume que $X_{k+1}$ es independiente de $S_k,$ las probabilidades de estos dos eventos componentes se calculan mediante multiplicación para obtener

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

Esta es la relación reportada por Stephan Kolassa en otro lugar en este hilo. Es una recursión de primer orden que comienza con $f_0 = 1$ (una suma vacía se define como cero, que ciertamente es par).

Hay muchas formas de resolverlo. Una de las más simples es escribir $q_k = 1 - 2p_k$ y expresar la recursión en términos de $a_k = 2f_k - 1$ como

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

Esta serie tiene la solución única $a_k = (1)(q_1)(q_2)\cdots(q_k),$ de donde (volviendo a la notación original en términos de $f_k$ y $p_k$)

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

exactamente la fórmula bellamente obtenida por Damian Pavlyshyn en otro lugar en este hilo.

Answer 3

Deberías poder calcular esto rápidamente a través de un enfoque de programación dinámica muy simple.

Sea $q_i$ la probabilidad de que $S_1+\dots+S_i$ sea par.

Entonces $q_1=1-p_1$.

Al pasar de $i-1$ a $i$, o bien tienes $S_1+\dots+S_{i-1}$ impar con probabilidad $1-q_{i-1}$, entonces $S_1+\dots+S_i$ será par si $S_i=1$, para una probabilidad total de $(1-q_{i-1})p_i$, o viceversa. En general,

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

Solo calcula $q_N$ iterando sobre $i$, y ahí lo tienes.

Answer 4

Aquí hay una solución al problema generalizado usando otro módulo

Las otras respuestas aquí son maravillosas, pero no puedo evitar generalizar tu problema. Considera el problema generalizado donde tenemos una secuencia $X_1,X_2,X_3,...$ con valores binarios independientes que tienen $X_i \sim \text{Bern}(p_i)$. Ahora define los residuos modulares de las sumas de interés como:

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

para $n=1,2,3,...$ y $k=2,3,4,...$. Queremos calcular la distribución de probabilidad para esta estadística, con probabilidades denotadas por:

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

Tu pregunta busca la probabilidad $\psi_{n,2}(0)$, que ocurre en el caso $k=2$ donde el valor $M_{n,2}$ es un indicador de que la suma de los primeros $n$ valores en la secuencia es par.

---

Computación recursiva de la probabilidad de interés: En este caso general podemos usar una extensión del algoritmo recursivo dado en la maravillosa respuesta de Stephen Kolassa. Para hacer esto, podemos usar la ecuación base:

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

y la ecuación recursiva:

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

Puedes encontrar una generalización similar relacionada con un problema asociado en esta respuesta relacionada. Es posible obtener una solución en forma cerrada para la función de probabilidad de interés (en lugar de una fórmula recursiva) utilizando el método descrito en la pregunta vinculada. Aquí procederemos utilizando la computación recursiva, lo que permite un cálculo eficiente de matrices de probabilidades en casos donde estamos interesados en todos los tamaños de muestra hasta un límite superior.

---

Implementación: Podemos implementar este método recursivo en R usando una función que tome n, k y p como entradas y genere el vector de probabilidades sobre todos los valores relevantes de $m$. (Nota: En la función a continuación, para la entrada p podemos tomarlo como un vector de probabilidades que tiene al menos $n$ elementos, o como una función que genere una probabilidad para cualquier número entero positivo para permitir una secuencia completa de probabilidades). La función se ajusta automáticamente al caso binario que te interesa, pero permite que el usuario establezca diferentes valores de k.

modular.prob <- function(n, p, k=2, keep.all = FALSE) {

  #Verificar las entradas n y k
  if (!is.numeric(n))      stop('Error: La entrada n debe ser un número')
  if (!is.numeric(k))      stop('Error: La entrada k debe ser un número')
  if (length(n) != 1)      stop('Error: La entrada n debe ser un solo número')
  if (length(k) != 1)      stop('Error: La entrada k debe ser un solo número')
  if (n != as.integer(n))  stop('Error: La entrada n debe ser un entero')
  if (k != as.integer(k))  stop('Error: La entrada k debe ser un entero')
  if (min(n < 0))          stop('Error: La entrada n debe ser un entero no negativo')
  if (min(k < 1))          stop('Error: La entrada n debe ser un entero positivo')

  #Verificar la entrada p
  if (is.function(p)) {
    if (length(formals(p)) != 1) {
      stop('Error: Si la entrada p es una función, solo debería tener un argumento') }
    pp <- rep(0, n)
    for (i in 1:n) { pp[i] <- p(i) }
  } else {
    if (!is.numeric(p))    stop('Error: La entrada p debe ser un vector numérico')
    if (length(p) < n)     stop('Error: La entrada p debe tener al menos n elementos')
    pp <- p[1:n]
    if (min(pp) < 0)       stop('Error: La entrada p debe tener valores entre cero y uno')
    if (max(pp) > 1)       stop('Error: La entrada p debe tener valores entre cero y uno') }

  #Verificar la entrada keep.all
  if (!is.logical(keep.all)) stop('Error: La entrada keep.all debe ser un valor lógico')
  if (length(keep.all) != 1) stop('Error: La entrada keep.all debe ser un solo valor lógico')

  #Establecer la matriz de probabilidades
  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('Par', 'Impar') }

  #Calcular probabilidades
  #Caso trivial donde k = 1 se trata sin recursión
  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))] } }

   #Dar salida
   if (keep.all) { PROBS } else { PROBS[n+1, ] } }

Aquí tienes un ejemplo que implementa la función con un vector de probabilidades específico con $n=20$ elementos. Primero calculamos las probabilidades para sumas pares/impares usando $k=2$ y luego calculamos las probabilidades para los restos modulares con $k=3$. En cada caso obtenemos una matriz de probabilidades para cada prueba y cada resto modular.

#Establecer vector de probabilidades (con veinte elementos)
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)

#Calcular probabilidades de sumas pares/impares
modular.prob(n = 20, p = p, keep.all = TRUE)

           Par       Impar
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

#Calcular probabilidades para restos modulares ternarios
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

Como se muestra en las salidas anteriores, típicamente veremos que $M_{n,k}$ converge en distribución a la distribución uniforme conforme $n \rightarrow \infty$. Esto se puede demostrar formalmente utilizando el TCL - ocurre debido a que la suma estandarizada de los valores binarios converge en distribución a una distribución continua (la distribución normal) y por lo tanto, el residuo modular sobre intervalos cada vez más cortos converge a probabilidades uniformes.