Éxito de los ensayos de Bernoulli con diferentes probabilidades.

Si se realizan 20 ensayos independientes de Bernoulli, cada uno con una probabilidad diferente de éxito y, por lo tanto, de fracaso. ¿Cuál es la probabilidad de que exactamente n de los 20 ensayos haya tenido éxito?

¿Existe una mejor manera de calcular estas probabilidades en lugar de simplemente sumar las combinaciones de probabilidades de éxito y fracaso?

La distribución por la que está preguntando se llama distribución binomial de Poisson , con un PMF bastante complicado (consulte Wikipedia para obtener una descripción más amplia)

$$\Pr(X=x) = \sum\limits_{A\in F_x} \prod\limits_{i\in A} p_i \prod\limits_{j\in A^c} (1-p_j)$$

Generalmente, el problema es que no puede usar esta ecuación para un número mayor de ensayos (generalmente cuando el número de ensayos excede ). También hay otros métodos para calcular el pmf, por ejemplo, fórmulas recursivas, pero son numéricamente inestables. La forma más fácil de solucionar esos problemas son los métodos de aproximación (descritos, por ejemplo, por Hong, 2013 ). Si definimos$n=30$

$$\mu = \sum_{i=1}^n p_i$$

$$\sigma = \sqrt{ \sum_{i=1}^n p_i(1-p_i) }$$

$$\gamma = \sigma^{-3} \sum_{i=1}^n p_i (1 - p_i) (1 - 2p_i)$$

entonces podemos aproximar pmf con distribución de Poisson a través de la ley de números pequeños o el teorema de Le Cams

$$\Pr(X = x) \approx \frac{\mu^x \exp(-\mu)}{x!}$$

pero ve que, en general, la aproximación binomial se comporta mejor ( Choi y Xia, 2002 )

$$\Pr(X = x) \approx \mathrm{Binom} \left( n, \frac{\mu}{n} \right)$$

puedes usar la aproximación normal

$$f(x) \approx \phi \left( \frac{x + 0.5 - \mu}{ \sigma} \right)$$

o cdf se puede aproximar utilizando la llamada aproximación normal refinada (Volkova, 1996)

$$F(x) \approx \max\left(0, \ g \left( \frac{x + 0.5 - \mu}{ \sigma} \right) \right)$$

donde .$g(x) = \Phi(x) + \gamma(1-x^2) \frac{\phi(x)}{6}$

Otra alternativa es, por supuesto, una simulación de Monte Carlo.

La dpbinomfunción R simple sería

dpbinom <- function(x, prob, log = FALSE,
                    method = c("MC", "PA", "NA", "BA"),
                    nsim = 1e4) {

  stopifnot(all(prob >= 0 & prob <= 1))
  method <- match.arg(method)

  if (method == "PA") {
    # poisson
    dpois(x, sum(prob), log)
  } else if (method == "NA") {
    # normal
    dnorm(x, sum(prob), sqrt(sum(prob*(1-prob))), log)
  } else if (method == "BA") {
    # binomial
    dbinom(x, length(prob), mean(prob), log)
  } else {
    # monte carlo
    tmp <- table(colSums(replicate(nsim, rbinom(length(prob), 1, prob))))
    tmp <- tmp/sum(tmp)
    p <- as.numeric(tmp[as.character(x)])
    p[is.na(p)] <- 0

    if (log) log(p)
    else p 
  }
}

La mayoría de los métodos (y más) también se implementan en el paquete R poibin .

---

Chen, LHY (1974). Sobre la convergencia de las distribuciones de Poisson Binomial a Poisson. Los Anales de Probabilidad, 2 (1), 178-180.

Chen, SX y Liu, JS (1997). Aplicaciones estadísticas de las distribuciones de Poisson-Binomial y condicional de Bernoulli. Statistica Sinica 7, 875-892.

Chen, SX (1993). Distribución de Poisson-Binomial, distribución condicional de Bernoulli y entropía máxima. Reporte técnico. Departamento de Estadística, Universidad de Harvard.

Chen, XH, Dempster, AP y Liu, JS (1994). Muestreo de población finita ponderada para maximizar la entropía. Biometrika 81, 457-469.

Wang, YH (1993). Sobre el número de éxitos en ensayos independientes. Statistica Sinica 3 (2): 295-312.

Hong, Y. (2013). Al calcular la función de distribución para la distribución binomial de Poisson. Estadística computacional y análisis de datos, 59, 41-51.

Volkova, AY (1996). Un refinamiento del teorema del límite central para sumas de indicadores aleatorios independientes. Teoría de la probabilidad y sus aplicaciones 40, 791-794.

Choi, KP y Xia, A. (2002). Aproximación del número de éxitos en ensayos independientes: Binomial versus Poisson. Los Anales de Probabilidad Aplicada, 14 (4), 1139-1148.

Le Cam, L. (1960). Un teorema de aproximación para la distribución binomial de Poisson. Pacific Journal of Mathematics 10 (4), 1181–1197.

Un enfoque es utilizar funciones generadoras. La solución a su problema es el coeficiente en el polinomio$x^n$

$$\prod_{i=1}^{20}(p_ix + 1-p_i).$$

Este es el equivalente de programación dinámica (tiempo cuadrático en el número de variables de Bernoulli) de hacer la suma en la distribución binomial de Poisson a partir de la respuesta de Tim (que sería el tiempo exponencial).