Como posso modelar eficientemente a soma das variáveis ​​aleatórias de Bernoulli?

38

YXipiY=XiPr(Xi=1)=piPr(Xi=0)=1pi .

Estou interessado em responder rapidamente a perguntas como (ondePr(Y<=k)k é fornecido).

Atualmente, uso simulações aleatórias para responder a essas perguntas. Eu desenho aleatoriamente cada Xi acordo com o seu pi , depois soma todos os valores de Xi para obter Y . Repito esse processo alguns milhares de vezes e retorno a fração de vezes Pr(Yk) .

Obviamente, isso não é totalmente preciso (embora a precisão aumente bastante à medida que o número de simulações aumenta). Além disso, parece que tenho dados suficientes sobre a distribuição para evitar as simulações de uso. Você consegue pensar em uma maneira razoável de obter a probabilidade exata Pr(Yk) ?

ps

Eu uso Perl & R.

EDITAR

Seguindo as respostas, pensei que poderiam ser necessários alguns esclarecimentos. Descreverei em breve a configuração do meu problema. Dado é um genoma circular com circunferência ce um conjunto de nfaixas mapeadas para ele. Por exemplo, c=3*10^9eranges={[100,200],[50,1000],[3*10^9-1,1000],...} . Observe que todos os intervalos estão fechados (ambas as extremidades são inclusivas). Observe também que lidamos apenas com números inteiros (unidades inteiras).

Estou procurando por regiões no círculo que sejam disfarçadas pelos nintervalos mapeados fornecidos . Portanto, para testar se um determinado intervalo de comprimento xno círculo é disfarçado, testo a hipótese de que os nintervalos são mapeados aleatoriamente. A probabilidade de um intervalo de comprimento mapeado q>xcobrir completamente o intervalo de comprimento especificado xé (q-x)/c. Essa probabilidade se torna bastante pequena quando cé grande e / ou qé pequena. O que me interessa é o número de intervalos (fora n) abrangidos x. É assimY se forma.

Testo minha hipótese nula vs. alternativa unilateral (disfarce). Observe também que estou testando várias hipóteses ( xcomprimentos diferentes ) e corrija isso.

David B
fonte
Seu p_i foi corrigido durante o exercício de modelagem ou eles podem mudar de um cálculo para o outro?
whuber
Os p_is são fixos.
David B
À luz das respostas atuais, você poderia compartilhar estimativas de (a) a soma dos p e (b) a soma de seus quadrados? Esses valores determinam suas opções.
whuber
@ Whuber: estes variam muito entre os casos. Não é um módulo único que estou criando (infelizmente).
David B
@ David Mas você não pode dar algumas orientações, como intervalos típicos? Por exemplo, se a soma dos p's varia entre 1 e 100, são informações úteis e sugerem algumas soluções eficientes, mas se podem chegar a 10.000, isso pode excluir algumas abordagens.
whuber

Respostas:

24

Se muitas vezes se assemelha a um Poisson , você tentou aproximar-se de um Poisson com o parâmetro λ=pi ?

EDIT : Encontrei um resultado teórico para justificar isso, bem como um nome para a distribuição de : é chamada de distribuição binomial de Poisson . A desigualdade de Le Cam mostra o quão próxima sua distribuição é aproximada pela distribuição de um Poisson com o parâmetro λ = p i . Diz-lhe a qualidade deste aprox é regido pela soma dos quadrados do p i é, parafraseando Steele (1994) . Então, se todos os seus p i s são razoavelmente pequeno, como agora parece que eles são, deve ser uma boa aproximação bonita.Yλ=pipipi

EDIT 2 : Quão pequeno é 'razoavelmente pequeno'? Bem, isso depende de quão boa você precisa da aproximação! O artigo da Wikipedia sobre o teorema de Le Cam fornece a forma precisa do resultado que mencionei acima: a soma das diferenças absolutas entre a função de massa de probabilidade (pmf) de e o pmf da distribuição de Poisson acima não é mais do que o dobro da soma das praças da p i s. Outro resultado de Le Cam (1960) pode ser mais fácil de usar: essa soma também não é superior a 18 vezes a maior p i . Existem mais alguns desses resultados ... ver Serfling (1978) para uma revisão.Ypipi

uma parada
fonte
1
+1 Não é uma má ideia. É provável que uma pequena mistura de Poissons faça um bom trabalho, dependendo de como a pergunta seja esclarecida.
whuber
1
Pensei em sugerir uma distribuição binomial negativa, que surge como uma mistura Gamma-Poisson, mas que apresenta uma variação maior que sua média, enquanto esse problema tem uma variação menor que sua média. Com base nisso, não tenho certeza se alguma mistura de Poissons funcionará, pois certamente essa mistura terá uma variação maior que a média?
onestop
@onestop Onde foi dito que a variação é menor que a média? Eu perdi essa afirmação.
whuber
Desculpe whuber, isso foi um pouco enigmático, mas esses comentários não permitem muita elaboração. das mpiktas é a variância, o que é menos do que a média, Σ p i . Apenas um pouco menos se o p i 's são em média muito pequena, porém, assim padrão Poisson pode ser uma boa o suficiente aprox. Talvez eu deva expandir minha resposta acima .. mas, em seguida, a conversa fica confusa. Bn=pi(1pi)pipi
onestop
O que você quer dizer com ? Como obtenho valores X i ? XiXi
David B
11

Encontrei sua pergunta enquanto procurava uma solução para esse mesmo problema. Não fiquei terrivelmente satisfeito com as respostas aqui, mas acho que há uma solução bastante simples que fornece a distribuição exata e é bastante tratável.

A distribuição da soma de duas variáveis ​​aleatórias discretas é a convolução de suas densidades. Portanto, se você tem onde conhece P ( X ) e P ( Y ) , pode calcular:Z=X+YP(X)P(Y)

P(Z=z)=k=P(X=k)P(Y=zk)

(Obviamente, para variáveis ​​aleatórias de Bernoulli, você não precisa ir muito ao infinito.)

Você pode usar isso para encontrar a distribuição exata da soma de seus RVs. Primeiro, some dois dos RVs, reunindo seus PDFs (por exemplo, [0,3, 0,7] * [0,6, 0,4] = [0,18, 0,54, 0,28]). Em seguida, convolva essa nova distribuição com seu próximo PDF de Bernoulli (por exemplo, [0,18, 0,54, 0,28] * [0,5, 0,5] = [0,09, 0,36, 0,41, 0,14]). Continue repetindo isso até que todos os RVs tenham sido adicionados. E pronto, o vetor resultante é o PDF exato da soma de todas as suas variáveis.

Eu verifiquei com simulação que isso produz os resultados corretos. Ele não se baseia em nenhuma suposição assintótica e não exige que os probes de Bernoulli sejam pequenos.

Também pode haver alguma maneira de fazer isso com mais eficiência do que convoluções repetidas, mas não pensei muito nisso. Espero que isso seja útil para alguém!

alex
fonte
2
Você já tentou isso com 40K variáveis ​​?? (Gostaria de saber quantas horas ou dias de cálculo são necessários ...)
whuber
5
(+1) Encontrei uma maneira de fazer essa ideia funcionar. Requer duas técnicas: primeiro, use FFT para as convoluções; segundo, não faça-os sequencialmente, mas divida e conquiste: faça-os em pares separados, depois faça os resultados em pares separados, etc. O algoritmo agora escala como vez de O ( n 2 ) para n probabilidades. Por exemplo, o Mathematica pode calcular toda a distribuição para 40.000 probabilidades em apenas 0,4 segundos. (1.000.000 são calculados em 10,5 segundos.) Fornecerei o código em um comentário subsequente. O(nlogn)O(n2)n
whuber
7
Aqui está o código do Mathematica : multinomial[p_] := Module[{lc, condense}, lc = Function[{s}, ListConvolve[s[[1]], s[[2]], {1, -1}, 0]]; condense = Function[{s}, Map[lc, Partition[s, 2, 2, {1, 1}, {{1}}]]]; Flatten[NestWhile[condense, Transpose[{1 - p, p}], Length[#] > 1 &]]] Para aplicá-lo, faça algo parecido p = RandomReal[{0, 1}, 40000]; pp = multinomial[p];. Isso cria as probabilidades pe calcula a distribuição exata pp. NB Quando a média de pnão é extrema, a distribuição está muito próxima do normal: isso leva a um algoritmo muito mais rápido ainda.
whuber
9

O @onestop fornece boas referências. O artigo da Wikipedia sobre distribuição binomial de Poisson fornece uma fórmula recursiva para calcular a distribuição exata de probabilidade; requer esforço de 2 ) . Infelizmente, é uma soma alternada, portanto será numericamente instável: não há como fazer isso com aritmética de ponto flutuante. Felizmente, quando o p i são pequenos, você só precisa calcular um pequeno número de probabilidades, de modo que o esforço é realmente proporcional ao O ( n log ( Σ i p i ) ) . A precisão necessária para realizar o cálculo com aritmética racional (O(n2)piO(nlog(ipi))isto é, exatamente, . Isso é viável.para que a instabilidade numérica não seja um problema) cresça lentamente o suficiente para que o tempo geral ainda seja aproximadamente O(n2)

Como teste, criei uma matriz de probabilidades para vários valores de n até n = 2 16 , que é o tamanho desse problema. Para pequenos valores de n (até n = 2 12 ), o tempo para o cálculo exato das probabilidades era em segundos e escalado quadraticamente, então arrisquei um cálculo para n = 2 16pi=1/(i+1)nn=216nn=212n=216até três DPs acima da média (probabilidades de 0, 1, ..., 22 sucessos). Demorou 80 minutos (com o Mathematica 8), de acordo com o tempo previsto. (As probabilidades resultantes são frações cujos numeradores e denominadores têm cerca de 75.000 dígitos cada um!) Isso mostra que o cálculo pode ser feito.

Uma alternativa é executar uma simulação longa (um milhão de tentativas deve ser feito). Ele só tem que ser feito uma vez, porque o não mudam.pi

whuber
fonte
9

(Como essa abordagem é independente das outras soluções postadas, incluindo uma que eu publiquei, estou oferecendo-a como uma resposta separada).

Você pode calcular a distribuição exata em segundos (ou menos), desde que a soma dos p's seja pequena.

Já vimos sugestões de que a distribuição possa ser aproximadamente gaussiana (em alguns cenários) ou Poisson (em outros cenários). De qualquer maneira, sabemos que sua média é a soma de p i e sua variância σ 2 é a soma de p i ( 1 - p i ) . Portanto, a distribuição será concentrada dentro de alguns desvios padrão de sua média, digamos z SDs com z entre 4 e 6 ou aproximadamente. Portanto, precisamos calcular apenas a probabilidade de que a soma X seja igual (um número inteiro) k para k = μμpiσ2pi(1pi)zzXk através de k = μ + z σ . Quando a maior parte do p ik=μzσk=μ+zσpi são pequeno, é aproximadamente igual a (mas ligeiramente menor do que) μ , de modo a ser conservadora que pode fazer o cálculo para k no intervalo [ μ - z σ2μk. Por exemplo, quando a soma depié igual a9e a escolha dez=6para cobrir bem as caudas, precisaríamos do cálculo para cobrirkem[9-6[μzμ,μ+zμ]pi9z=6k=[0,27], que são apenas 28 valores.[969,9+69][0,27]

A distribuição é calculada recursivamente . Vamos ser a distribuição da soma do primeiro i dessas variáveis Bernoulli. Para qualquer j de 0 a i + 1 , a soma das primeiras variáveis i + 1 pode ser igual a j de duas maneiras mutuamente exclusivas: a soma das primeiras variáveis i é igual a j e o i + 1 st é 0 ou então a soma de a primeira i variáveis ​​é igual a j - 1 e afiij0i+1i+1jiji+1st0ij1i+1st is 1. Therefore

fi+1(j)=fi(j)(1pi+1)+fi(j1)pi+1.

We only need to carry out this computation for integral j in the interval from max(0,μzμ) to μ+zμ.

When most of the pi are tiny (but the 1pi are still distinguishable from 1 with reasonable precision), this approach is not plagued with the huge accumulation of floating point roundoff errors used in the solution I previously posted. Therefore, extended-precision computation is not required. For example, a double-precision calculation for an array of 216 probabilities pi=1/(i+1) (μ=10.6676, requiring calculations for probabilities of sums between 0 and 31) took 0.1 seconds with Mathematica 8 and 1-2 seconds with Excel 2002 (both obtained the same answers). Repeating it with quadruple precision (in Mathematica) took about 2 seconds but did not change any answer by more than 3×1015. Terminating the distribution at z=6 SDs into the upper tail lost only 3.6×108 of the total probability.

Another calculation for an array of 40,000 double precision random values between 0 and 0.001 (μ=19.9093) took 0.08 seconds with Mathematica.

This algorithm is parallelizable. Just break the set of pi into disjoint subsets of approximately equal size, one per processor. Compute the distribution for each subset, then convolve the results (using FFT if you like, although this speedup is probably unnecessary) to obtain the full answer. This makes it practical to use even when μ gets large, when you need to look far out into the tails (z large), and/or n is large.

The timing for an array of n variables with m processors scales as O(n(μ+zμ)/m). Mathematica's speed is on the order of a million per second. For example, with m=1 processor, n=20000 variates, a total probability of μ=100, and going out to z=6 standard deviations into the upper tail, n(μ+zμ)/m=3.2 million: figure a couple seconds of computing time. If you compile this you might speed up the performance two orders of magnitude.

Incidentally, in these test cases, graphs of the distribution clearly showed some positive skewness: they aren't normal.

For the record, here is a Mathematica solution:

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 The color coding applied by this site is meaningless for Mathematica code. In particular, the gray stuff is not comments: it's where all the work is done!)

An example of its use is

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

Edit

An R solution is ten times slower than Mathematica in this test case--perhaps I have not coded it optimally--but it still executes quickly (about one second):

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)

Plot of PDF

whuber
fonte
8

With different pi your best bet I think is normal approximation. Let Bn=i=1npi(1pi). Then

Bn1/2(i=1nXii=1npi)N(0,1),
as n, provided that for each ε>0

Bn1i=1nE((Xipi)21{|Xipi|>εBn1/2})0,
as n, which for Bernoulli variables will hold if Bn. This is the so-called Lindeberg condition, which is sufficient and necessary for convergence to the standard normal.

Update: The approximation error can be calculated from the following inequality:

supx|Fn(x)Φ(x)|ALn,
where
Ln=Bn3/2i=1nE|Xipi|3
and Fn is the cdf of the scaled and centered sum of Xi.

As whuber pointed out, the convergence can be slow for badly behaved pi. For pi=11+i we have Bnlnn and Ln(lnn)1/2. Then taking n=216 we get that the maximum deviation from the standard normal cdf is a whopping 0.3.

mpiktas
fonte
3
This is not true when the p_i approach zero as i increases. Otherwise, you have just proven that the Poisson distribution is Normal!
whuber
1
That is why it must be Bn. If pi approach zero at rate faster than 1/i, limBn<.
M101
@mpiktas is right. The analogy to the Poisson distribution doesn't quite fit, here.
By the way, I didn't actually check that monstrous condition in the second paragraph.
@G. Jay Kerns I agree that the analogy to the Poisson is imperfect, but I think it gives good guidance. Imagine a sequence of p's, p_i = 10^{-j}, where j is the order of magnitude of i (equal to 1 for i <= 10, to 2 for i <= 100, etc.). When n = 10^k, 90% of the p's equal 10^{-k} and their sum looks Poisson with expectation 0.9. Another 9% equal 10^{1-k} and their sum looks Poisson (with the same expectation). Thus the distribution looks approximately like a sum of k Poisson variates. It's obviously nowhere near Normal. Whence the need for the "monstrous condition."
whuber
4

Well, based on your description and the discussion in the comments it is clear that Y has mean ipi and variance ipi(1pi). The shape of Y's distribution will ultimately depend on the behavior of pi. For suitably "nice" pi (in the sense that not too many of them are really close to zero), the distribution of Y will be approximately normal (centered right at pi). But as ipi starts heading toward zero the distribution will be shifted to the left and when it crowds up against the y-axis it will start looking a lot less normal and a lot more Poisson, as @whuber and @onestop have mentioned.

From your comment "the distribution looks Poisson" I suspect that this latter case is what's happening, but can't really be sure without some sort of visual display or summary statistics about the p's. Note however, as @whuber did, that with sufficiently pathological behavior of the p's you can have all sorts of spooky things happen, like limits that are mixture distributions. I doubt that is the case here, but again, it really depends on what your p's are doing.

As to the original question of "how to efficiently model", I was going to suggest a hierarchical model for you but it isn't really appropriate if the p's are fixed constants. In short, take a look at a histogram of the p's and make a first guess based on what you see. I would recommend the answer by @mpiktas (and by extension @csgillespie) if your p's aren't too crowded to the left, and I would recommend the answer by @onestop if they are crowded left-ly.

By the way, here is the R code I used while playing around with this problem: the code isn't really appropriate if your p's are too small, but it should be easy to plug in different models for p (including spooky-crazy ones) to see what happens to the ultimate distribution of 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)))

Now take a look at the results.

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

Have fun; I sure did.


fonte
Why do you say "the code isn't really appropriate if your ps are too small"? Seems to work ok to me, e.g. with shape1=1, shape2=999, giving a mean p of 0.001.
onestop
@onestop what I meant was the specific choice of (1,10) written above doesn't give values of p that are very small, to the point that the normal approximation looks pretty good. If a person wanted the Poisson to come out then they would need to try something else; it sounds like your choice of (1,999) does a good job, yes? I had also thought to make α<1, say, 0.25, but I haven't tried that.
2

I think other answers are great, but I didn't see any Bayesian ways of estimating your probability. The answer doesn't have an explicit form, but the probability can be simulated using R.

Here is the attempt:

Xi|piBer(pi)

piBeta(α,β)

Using wikipedia we can get estimates of α^ and β^ (see parameter estimation section).

Now you can generate draws for the ith step, generate pi from Beta(α^,β^) and then generate Xi from Ber(pi). After you have done this N times you can get Y=Xi. This is a single cycle for generation of Y, do this M(large) number of times and the histogram for M Ys will be the estimate of density of Y.

Prob[Yy]=#YyM

This analysis is valid only when pi are not fixed. This is not the case here. But I will leave it here, in case someone has a similar question.

suncoolsu
fonte
1
To some purists this may not be Bayesian. This is actually empirical Bayesian, but it is a quick way to simulate your probabilities in R, without resorting to hyper prior mumbo jumbo.
suncoolsu
1
Why do you need priors when the p_i are given?
whuber
@whuber. Thanks, you are right. I missed the fixed part. I thought David is just using the value to be pi as (q-x)/c and is not fixed. I will edit my answer.
suncoolsu
@suncoolsu - note that a "beta-bernoulli" distribution is just another bernoulli distribution but replacing piαα+β. This is becase (1xi)B(α+xi,β+1xi)B(α,β)=αxiβ1xiα+β. So basically by mixing over pi you are applying the binomial approximation here p1=p2==pn.
probabilityislogic
2

As has been mentioned in other answers, the probability distribution you describe is the Poisson Binomial distribution. An efficient method for computing the CDF is given in Hong, Yili. On computing the distribution function for the Poisson binomial distribution.

The approach is to efficiently compute the DFT (discrete Fourier transform) of the characteristic function.

The characteristic function of the Poisson binomial distribution is give by ϕ(t)=jn[(1pj)+pjeit] (i=1).

The algorithm is:

  1. Let zj(k)=1pj+pjcos(ωk)+ipjsin(ωk), for ω=2πn+1.
  2. Define xk=exp{jnlog(zj(k))}, define x0=1.
  3. Compute xk for k=1,,[n/2]. Use symmetry x¯k=xn+1k to get the rest.
  4. Apply FFT to the vector 1n+1<x0,x1,,xn>.
  5. Take the cumulative sum of result to get the CDF.

The algorithm is available in the poibin R package.

This approach gives much better results than the recursive formulations as they tend to lack numerical stability.

Kyle
fonte
3
I have access only to the abstract of that paper, but it sounds like it implements the method I used at stats.stackexchange.com/questions/41247/… and discusses how it performs compares to the other methods given in this thread. If you know more about what the paper has accomplished, we would be glad to read a summary.
whuber
1

I would suggest applying Poisson approximation. It is well known (see A. D. Barbour, L. Holst and S. Janson: Poisson Approximation) that the total variation distance between Y and a r.v. Z having Poisson distribution with the parameter ipi is small:

supA|P(YA)P(ZA)|min{1,1ipi}ipi2.
There are also bounds in terms of information divergence (the Kullback-Leibler distance, you may see P. Harremoёs: Convergence to the Poisson Distribution in Information Divergence. Preprint no. 2, Feb. 2003, Mathematical Department, University of Copenhagen. http://www.harremoes.dk/Peter/poisprep.pdf and other publications of P.Harremoёs), chi-squared distance (see Borisov and Vorozheikin https://link.springer.com/article/10.1007%2Fs11202-008-0002-3) and some other distances.

For the accuracy of approximation |Ef(Y)Ef(Z)| for unbounded functions f you may see Borisov and Ruzankin https://projecteuclid.org/euclid.aop/1039548369 . Besides, that paper contains a simple bound for probabilities: for all A, we have

P(YA)1(1maxipi)2P(ZA).

Pavel Ruzankin
fonte
1
+1 Thank you for the useful quantitative information about the approximation bounds. Welcome to our site!
whuber