¿Cómo puedo modelar eficazmente la suma de variables aleatorias Bernoulli?

Question

Estoy modelando una variable aleatoria ( $Y$ ) que es la suma de unas ~15-40k variables aleatorias Bernoulli independientes ( $X_i$ ), cada una con una probabilidad de éxito diferente ( $p_i$ ). Formalmente, $Y=\sum X_i$ donde $\Pr(X_i=1)=p_i$ y $\Pr(X_i=0)=1-p_i$ .

Estoy interesado en responder rápidamente a preguntas como $\Pr(Y<=k)$ (donde $k$ se da).

Actualmente, utilizo simulaciones aleatorias para responder a estas preguntas. Saco al azar cada $X_i$ según su $p_i$ , entonces suma todos los $X_i$ para obtener $Y'$ . Repito este proceso unos cuantos miles de veces y devuelvo la fracción de veces $\Pr(Y'\leq k)$ .

Obviamente, esto no es totalmente preciso (aunque la precisión aumenta en gran medida a medida que aumenta el número de simulaciones). Además, parece que tengo suficientes datos sobre la distribución para evitar las simulaciones de uso. ¿Se te ocurre alguna forma razonable de obtener la probabilidad exacta $\Pr(Y\leq k)$ ?

p.d.

Utilizo Perl y R.

EDITAR

A raíz de las respuestas, he pensado que sería necesario hacer algunas aclaraciones. En breve describiré el escenario de mi problema. Dado es un genoma circular con circunferencia c y un conjunto de n rangos asignados a ella. Por ejemplo, c=3*10^9 y ranges={[100,200],[50,1000],[3*10^9-1,1000],...} . Observe que todos los rangos son cerrados (ambos extremos son inclusivos). También hay que tener en cuenta que sólo tratamos con números enteros (unidades enteras).

Estoy buscando regiones en el círculo que no están cubiertas por el n rangos mapeados. Así que para comprobar si un rango dado de longitud x en el círculo está subcubierto, pruebo la hipótesis de que el n Los rangos se asignan de forma aleatoria. La probabilidad de que un rango mapeado de longitud q>x cubrirá completamente el rango de longitud dado x es (q-x)/c . Esta probabilidad se hace bastante pequeña cuando c es grande y/o q es pequeño. Lo que me interesa es el número de rangos (de n ) que cubren x . Así es como Y se forma.

Pongo a prueba mi hipótesis nula frente a la alternativa unilateral (infracobertura). También hay que tener en cuenta que estoy probando múltiples hipótesis (diferentes x longitudes), y asegúrese de corregirlo.

Answer 1

Si a menudo se parece a un Poisson ¿has probado a aproximarla por una Poisson con parámetro $\lambda = \sum p_i$ ?

EDITAR : He encontrado un resultado teórico que lo justifica, así como un nombre para la distribución de $Y$ : se llama el Distribución binomial de Poisson . La desigualdad de Le Cam le indica el grado de aproximación de su distribución a la distribución de una Poisson con parámetro $\lambda = \sum p_i$ . Te dice que la calidad de este aprox se rige por la suma de los cuadrados de los $p_i$ s, parafraseando Steele (1994) . Así que si todo su $p_i$ s son razonablemente pequeños, como parece ser ahora, debería ser una aproximación bastante buena.

EDITAR 2 : ¿Cómo de pequeño es "razonablemente pequeño"? Bueno, eso depende de lo buena que sea la aproximación. El Artículo de Wikipedia sobre el teorema de Le Cam da la forma precisa del resultado al que me refería antes: la suma de las diferencias absolutas entre los función de masa de probabilidad (pmf) de $Y$ y la pmf de la distribución de Poisson anterior no es más que el doble de la suma de los cuadrados de las $p_i$ s. Otro resultado de Le Cam (1960) puede ser más fácil de usar: esta suma tampoco es más de 18 veces la mayor $p_i$ . Hay bastantes más resultados de este tipo... ver Serfling (1978) para una revisión.

Answer 2

Me encontré con su pregunta mientras buscaba una solución a este mismo problema. No quedé muy satisfecho con las respuestas aquí, pero creo que hay una solución bastante sencilla que te da la distribución exacta, y es bastante manejable.

La distribución de la suma de dos variables aleatorias discretas es la convolución de sus densidades. Así, si se tiene $Z = X + Y$ donde se sabe $P(X)$ y $P(Y)$ entonces puedes calcular:

$$P(Z=z) = \sum_{k=-\infty}^{\infty} P(X=k) \; P(Y=z-k)$$

(Por supuesto, para las variables aleatorias de Bernoulli no es necesario ir bastante hasta el infinito).

Puedes usar esto para encontrar la distribución exacta de la suma de tus VR. En primer lugar, sume dos de los VR mediante la convolución de sus PDF (por ejemplo, [0,3, 0,7] * [0,6, 0,4] = [0,18, 0,54, 0,28]). A continuación, convierta esa nueva distribución con su siguiente PDF Bernoulli (por ejemplo, [0,18, 0,54, 0,28] * [0,5, 0,5] = [0,09, 0,36, 0,41, 0,14]). Repite esto hasta que se hayan añadido todos los VR. Y voilá, el vector resultante es la PDF exacta de la suma de todas sus variables.

He comprobado con la simulación que esto produce los resultados correctos. No se basa en ningún supuesto asintótico, y no tiene requisitos de que las probabilidades de Bernoulli sean pequeñas.

También puede haber alguna manera de hacer esto de manera más eficiente que la convolución repetida, pero no he pensado en ello muy profundamente. ¡Espero que esto sea útil para alguien!

Answer 3

(Dado que este enfoque es independiente de las otras soluciones publicadas, incluida una que yo he publicado, lo ofrezco como una respuesta separada).

Se puede calcular la distribución exacta en segundos (o menos) siempre que la suma de las p sea pequeña.

Ya hemos visto sugerencias de que la distribución podría ser aproximadamente gaussiana (en algunos escenarios) o de Poisson (en otros escenarios). En cualquier caso, sabemos que su media $\mu$ es la suma de los $p_i$ y su varianza $\sigma^2$ es la suma de $p_i(1-p_i)$ . Por lo tanto, la distribución se concentrará dentro de unas pocas desviaciones estándar de su media, digamos $z$ SDs con $z$ entre 4 y 6 o más. Por lo tanto, sólo tenemos que calcular la probabilidad de que la suma $X$ es igual a (un número entero) $k$ para $k = \mu - z \sigma$ a través de $k = \mu + z \sigma$ . Cuando la mayoría de los $p_i$ son pequeños, $\sigma^2$ es aproximadamente igual (pero ligeramente inferior) a $\mu$ por lo que, para ser conservadores, podemos hacer el cálculo para $k$ en el intervalo $[\mu - z \sqrt{\mu}, \mu + z \sqrt{\mu}]$ . Por ejemplo, cuando la suma de los $p_i$ es igual a $9$ y elegir $z = 6$ para cubrir bien las colas, necesitaríamos que el cómputo cubriera $k$ en $[9 - 6 \sqrt{9}, 9 + 6 \sqrt{9}]$ = $[0, 27]$ , que son sólo 28 valores.

La distribución se calcula recursivamente . Sea $f_i$ sea la distribución de la suma de los primeros $i$ de estas variables Bernoulli. Para cualquier $j$ de $0$ a través de $i+1$ la suma de la primera $i+1$ las variables pueden ser iguales a $j$ de dos maneras mutuamente excluyentes: la suma del primer $i$ variables es igual a $j$ y el $i+1^\text{st}$ es $0$ o bien la suma del primer $i$ variables es igual a $j-1$ y el $i+1^\text{st}$ es $1$ . Por lo tanto,

$$f_{i+1}(j) = f_i(j)(1 - p_{i+1}) + f_i(j-1) p_{i+1}.$$

Sólo tenemos que realizar este cálculo para la integral $j$ en el intervalo de $\max(0, \mu - z \sqrt{\mu})$ a $\mu + z \sqrt{\mu}.$

Cuando la mayoría de los $p_i$ son diminutos (pero el $1 - p_i$ siguen siendo distinguibles de $1$ con una precisión razonable), este enfoque no está plagado de la enorme acumulación de errores de redondeo en coma flotante utilizados en la solución que publiqué anteriormente. Por lo tanto, no es necesario el cálculo de precisión extendida. Por ejemplo, un cálculo de doble precisión para un array de $2^{16}$ probabilidades $p_i = 1/(i+1)$ ( $\mu = 10.6676$ , requiriendo cálculos para las probabilidades de las sumas entre $0$ y $31$ ) tardó 0,1 segundos con Mathematica 8 y 1-2 segundos con Excel 2002 (ambos obtuvieron las mismas respuestas). Repetirlo con precisión cuádruple (en Mathematica) llevó unos 2 segundos pero no cambió ninguna respuesta en más de $3 \times 10^{-15}$ . Terminar la distribución en $z = 6$ Los SD en la cola superior sólo perdieron $3.6 \times 10^{-8}$ de la probabilidad total.

Otro cálculo para una matriz de 40.000 valores aleatorios de doble precisión entre 0 y 0,001 ( $\mu = 19.9093$ ) tardó 0,08 segundos con Mathematica.

Este algoritmo es paralelizable. Sólo hay que romper el conjunto de $p_i$ en subconjuntos disjuntos de tamaño aproximadamente igual, uno por procesador. Calcule la distribución de cada subconjunto y, a continuación, convolucione los resultados (utilizando la FFT si lo desea, aunque esta aceleración es probablemente innecesaria) para obtener la respuesta completa. Esto hace que sea práctico utilizarlo incluso cuando $\mu$ se hace grande, cuando hay que mirar lejos en las colas ( $z$ grande), y/o $n$ es grande.

El tiempo para una serie de $n$ variables con $m$ escalas de procesadores como $O(n(\mu + z \sqrt{\mu})/m)$ . La velocidad de Mathematica es del orden de un millón por segundo. Por ejemplo, con $m = 1$ procesador, $n = 20000$ variantes, una probabilidad total de $\mu = 100$ y salir a $z = 6$ desviaciones estándar en la cola superior, $n(\mu + z \sqrt{\mu})/m = 3.2$ millones: calcula un par de segundos de tiempo de cálculo. Si lo compilas, podrías acelerar el rendimiento dos órdenes de magnitud.

Por cierto, en estos casos de prueba, los gráficos de la distribución mostraban claramente cierta asimetría positiva: no son normales.

Para que conste, aquí hay una solución de Mathematica:

pb[p_, z_] := Module[
  {\[Mu] = Total[p]},
  Fold[#1 - #2 Differences[Prepend[#1, 0]] &, 
   Prepend[ConstantArray[0, Ceiling[\[Mu] + Sqrt[\[Mu]] z]], 1], p]
  ]

( NB El código de colores aplicado por este sitio no tiene sentido para el código de Mathematica. En particular, lo gris es no comentarios: ¡es donde se hace todo el trabajo!)

Un ejemplo de su uso es

pb[RandomReal[{0, 0.001}, 40000], 8]

---

Editar

Un R es diez veces más lenta que Mathematica en este caso de prueba -quizás no lo he codificado de forma óptima- pero aún así se ejecuta rápidamente (alrededor de un segundo):

pb <- function(p, z) {
  mu <- sum(p)
  x <- c(1, rep(0, ceiling(mu + sqrt(mu) * z)))
  f <- function(v) {x <<- x - v * diff(c(0, x));}
  sapply(p, f); x  
}
y <- pb(runif(40000, 0, 0.001), 8)
plot(y)

Answer 4

@onestop proporciona buenas referencias. El artículo de Wikipedia sobre la distribución binomial de Poisson ofrece una fórmula recursiva para calcular la distribución de probabilidad exacta; requiere $O(n^2)$ esfuerzo. Por desgracia, se trata de una suma alternante, por lo que será numéricamente inestable: es inútil hacer este cálculo con aritmética de punto flotante. Afortunadamente, cuando el $p_i$ son pequeñas, sólo hay que calcular un pequeño número de probabilidades, por lo que el esfuerzo es realmente proporcional a $O(n \log(\sum_i{p_i}))$ . La precisión necesaria para realizar el cálculo con aritmética racional ( es decir, exactamente, para que la inestabilidad numérica no sea un problema) crece con la suficiente lentitud como para que el tiempo total pueda seguir siendo aproximadamente $O(n^2)$ . Eso es factible.

Como prueba, he creado una matriz de probabilidades $p_i = 1/(i+1)$ para varios valores de $n$ hasta $n = 2^{16}$ que es el tamaño de este problema. Para valores pequeños de $n$ (hasta $n = 2^{12}$ ) el tiempo para el cálculo exacto de las probabilidades estaba en segundos y se escalaba cuadráticamente, así que me aventuré a hacer un cálculo para $n = 2^{16}$ hasta tres DE por encima de la media (probabilidades de 0, 1, ..., 22 aciertos). Se tardó 80 minutos (con Mathematica 8), en línea con el tiempo previsto. (Las probabilidades resultantes son fracciones cuyos numeradores y denominadores tienen unos 75.000 dígitos cada uno). Esto demuestra que el cálculo puede hacerse.

Una alternativa es realizar una simulación larga (debería bastar con un millón de pruebas). Sólo hay que hacerlo una vez, porque el $p_i$ no cambian.

Answer 5

Con diferentes $p_i$ su mejor apuesta creo que es la aproximación normal. Dejemos que $B_n=\sum_{i=1}^np_i(1-p_i)$ . Entonces

\begin{align*} B_n^{-1/2}\left(\sum_{i=1}^nX_i-\sum_{i=1}^np_i\right)\to N(0,1), \end{align*} como $n\to\infty$ siempre que para cada $\varepsilon>0$

\begin{align*} B_n^{-1}\sum_{i=1}^nE\left((X_i-p_i)^2\mathbf{1}\{|X_i-p_i|>\varepsilon B_n^{1/2}\}\right)\to 0, \end{align*} como $n\to\infty$ que para las variables de Bernoulli se cumple si $B_n\to\infty$ . Esta es la llamada condición de Lindeberg, que es suficiente y necesaria para la convergencia a la normal estándar.

Actualización: El error de aproximación se puede calcular a partir de la siguiente desigualdad:

\begin{align*} \sup_x|F_n(x)-\Phi(x)|\le AL_n, \end{align*} donde \begin{align*} L_n=B_n^{-3/2}\sum_{i=1}^nE|X_i-p_i|^3 \end{align*} y $F_n$ es la fdc de la suma escalada y centrada de $X_i$ .

Como ha señalado Whuber, la convergencia puede ser lenta para los que se comportan mal $p_i$ . Para $p_i=\frac{1}{1+i}$ tenemos $B_n\approx \ln n$ y $L_n\approx (\ln n)^{-1/2}$ . Entonces, tomando $n=2^{16}$ obtenemos que la máxima desviación de la cdf normal estándar es la friolera de 0,3.