¿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

Pues bien, en base a su descripción y a la discusión en los comentarios está claro que $Y$ tiene media $\sum_i p_i$ y la varianza $\sum_i p_{i}(1-p_{i})$ . La forma de $Y$ dependerá en última instancia del comportamiento de $p_i$ . Para ser convenientemente "agradable" $p_i$ (en el sentido de que no hay demasiados cerca de cero), la distribución de $Y$ será aproximadamente normal (centrado justo en $\sum p_i$ ). Pero como $\sum_i p_i$ comienza a dirigirse hacia el cero la distribución se desplazará hacia la izquierda y cuando se agolpe contra el $y$ -eje comenzará a parecer mucho menos normal y mucho más Poisson, como han mencionado @whuber y @onestop.

Por tu comentario "la distribución parece de Poisson" sospecho que este último caso es el que está ocurriendo, pero no puedo estar seguro sin algún tipo de presentación visual o estadísticas resumidas sobre la $p$ 's. Nótese sin embargo, como hizo @whuber, que con un comportamiento suficientemente patológico del $p$ puedes tener todo tipo de cosas espeluznantes, como límites que son distribuciones de mezcla. Dudo que sea el caso aquí, pero de nuevo, realmente depende de lo que su $p$ de la UE.

En cuanto a la pregunta original de "cómo modelar eficientemente", iba a sugerirte un modelo jerárquico pero no es realmente apropiado si el $p$ son constantes fijas. En resumen, mira un histograma de la $p$ y hacer una primera conjetura basada en lo que ves. Yo recomendaría la respuesta de @mpiktas (y por extensión @csgillespie) si su $p$ No hay demasiada aglomeración a la izquierda, y yo recomendaría la respuesta de @onestop si están aglomerados a la izquierda.

Por cierto, aquí está el código R que utilicé mientras jugaba con este problema: el código no es realmente apropiado si su $p$ son demasiado pequeños, pero debería ser fácil conectar diferentes modelos para $p$ (incluidas las espeluznantes) para ver qué ocurre con la distribución final de $Y$ .

set.seed(1)
M <- 5000
N <- 15000
p <- rbeta(N, shape1 = 1, shape2 = 10)
Y <- replicate(M, sum(rbinom(N, size = 1, prob = p)))

Ahora echa un vistazo a los resultados.

hist(Y)
mean(Y)
sum(p)
var(Y)
sum(p*(1 - p))

Diviértete; yo lo hice.

Answer 2

Creo que otras respuestas son geniales, pero no he visto ninguna forma bayesiana de estimar su probabilidad. La respuesta no tiene una forma explícita, pero la probabilidad puede ser simulada usando R.

Aquí está el intento:

$$ X_i | p_i \sim Ber(p_i)$$

$$ p_i \sim Beta(\alpha, \beta) $$

Utilizando wikipedia podemos obtener estimaciones de $\hat{\alpha}$ y $\hat{\beta}$ (véase la sección de estimación de parámetros).

Ahora puede generar sorteos para el $i^{th}$ paso, generar $p_i$ de $Beta(\hat{\alpha},\hat{\beta})$ y luego generar $ X_i$ de $Ber(p_i)$ . Después de hacer esto $N$ veces se puede conseguir $Y = \sum X_i$ . Este es un ciclo único para la generación de Y, haga esto $M$ (gran) número de veces y el histograma para $M$ Ys será la estimación de la densidad de Y.

$$Prob[Y \leq y] = \frac {\#Y \leq y} {M}$$

Este análisis sólo es válido cuando $p_i$ no son fijos. Este no es el caso aquí. Pero lo dejaré aquí, por si alguien tiene una pregunta similar.

Answer 3

Yo sugeriría aplicar la aproximación de Poisson. Es bien sabido (véase A. D. Barbour, L. Holst y S. Janson: Poisson Approximation) que la distancia de variación total entre $Y$ y un r.v. $Z$ con una distribución de Poisson con el parámetro $\sum_i p_i$ es pequeño: $$ \sup_A |{\bf P}(Y\in A) - {\bf P}(Z\in A)| \le \min \left\{ 1, \frac{1}{\sum_i p_i} \right\} \sum_i p_i^2. $$ También hay límites en términos de divergencia de información (la distancia de Kullback-Leibler, puede ver P. Harremoёs: Convergencia a la distribución de Poisson en la divergencia de información. Preprint no. 2, Feb. 2003, Departamento de Matemáticas, Universidad de Copenhague. http://www.harremoes.dk/Peter/poisprep.pdf y otras publicaciones de P.Harremoёs), la distancia chi-cuadrado (véase Borisov y Vorozheikin https://link.springer.com/article/10.1007%2Fs11202-008-0002-3 ) y algunas otras distancias.

Para la precisión de la aproximación $|{\bf E}f(Y) - {\bf E}f(Z)|$ para funciones no limitadas $f$ puede ver a Borisov y Ruzankin https://projecteuclid.org/euclid.aop/1039548369 . Además, ese documento contiene un límite simple para las probabilidades: para todo $A$ tenemos $${\bf P}(Y\in A) \le \frac{1}{(1-\max_i p_i)^2} {\bf P}(Z\in A).$$

Answer 4

Como se ha mencionado en otras respuestas, la distribución de probabilidad que describes es la distribución Binomial de Poisson. Un método eficiente para calcular la FCD se da en Hong, Yili. Sobre el cálculo de la función de distribución para la distribución binomial de Poisson .

El enfoque consiste en calcular eficazmente la DFT (transformada discreta de Fourier) de la función característica.

La función característica de la distribución binomial de Poisson viene dada por $\phi(t) = \prod_j^n [(1-p_j)+p_je^{it}]$ ( $i=\sqrt{-1}$ ).

El algoritmo es:

1. Dejemos que $z_j(k) = 1-p_j+p_j \text{cos}(\omega k)+ i p_j \text{sin}(\omega k)$ , para $\omega=\frac{2\pi}{n+1}$ .
2. Definir $x_k=\text{exp}\{\sum_j^n log(z_j(k))\}$ , defina $x_0=1$ .
3. Calcula $x_k$ para $k=1,\dots,[n/2]$ . Utilizar la simetría $\bar{x}_k=x_{n+1-k}$ para conseguir el resto.
4. Aplicar la FFT al vector $\frac{1}{n+1}<x_0,x_1,\dots,x_n>$ .
5. Tome la suma acumulada del resultado para obtener la FCD.

El algoritmo está disponible en el poibin Paquete R.

Este enfoque da resultados mucho mejores que las formulaciones recursivas, ya que tienden a carecer de estabilidad numérica.