Bagaimana saya bisa memodelkan jumlah variabel acak Bernoulli secara efisien?

38

Saya memodelkan variabel acak ( ) yang merupakan jumlah dari beberapa ~ 15-40k variabel acak Bernoulli independen ( ), masing-masing dengan probabilitas keberhasilan yang berbeda ( ). Secara formal, mana dan \ Pr (X_i = 0) = 1-p_i .YXipiY=XiPr(Xi=1)=piPr(Xi=0)=1pi

Saya tertarik untuk dengan cepat menjawab pertanyaan seperti Pr(Y<=k) (di mana k diberikan).

Saat ini, saya menggunakan simulasi acak untuk menjawab pertanyaan seperti itu. Saya secara acak menggambar setiap Xi sesuai dengan p_i -nya pi, lalu menjumlahkan semua nilai Xi untuk mendapatkan Y . Saya ulangi proses ini beberapa ribu kali dan mengembalikan fraksi kali Pr(Yk) .

Jelas, ini tidak sepenuhnya akurat (walaupun keakuratannya meningkat dengan meningkatnya jumlah simulasi). Juga, sepertinya saya punya cukup data tentang distribusi untuk menghindari penggunaan simulasi. Bisakah Anda memikirkan cara yang masuk akal untuk mendapatkan probabilitas yang tepat Pr(Yk) ?

ps

Saya menggunakan Perl & R.

EDIT

Mengikuti tanggapan saya pikir beberapa klarifikasi mungkin diperlukan. Saya akan segera menjelaskan pengaturan masalah saya. Diberikan adalah genom melingkar dengan keliling cdan seperangkat nrentang yang dipetakan untuk itu. Sebagai contoh, c=3*10^9dan ranges={[100,200],[50,1000],[3*10^9-1,1000],...}. Perhatikan bahwa semua rentang ditutup (kedua ujungnya termasuk). Perhatikan juga bahwa kami hanya berurusan dengan bilangan bulat (seluruh unit).

Saya mencari daerah di lingkaran yang tertutup oleh nrentang yang dipetakan. Jadi untuk menguji apakah rentang panjang tertentu xpada lingkaran tertutup, saya menguji hipotesis bahwa nrentang tersebut dipetakan secara acak. Probabilitas rentang panjang yang dipetakan q>xakan sepenuhnya mencakup rentang panjang yang xdiberikan (q-x)/c. Probabilitas ini menjadi sangat kecil ketika cbesar dan / atau qkecil. Yang saya minati adalah jumlah rentang (dari n) yang mencakup x. Beginilah cara Yterbentuk.

Saya menguji hipotesis nol saya vs. alternatif satu sisi (undercoverage). Perhatikan juga saya sedang menguji beberapa hipotesis ( xpanjang yang berbeda ), dan pastikan untuk memperbaiki ini.

David B
sumber
Apakah p_i Anda diperbaiki selama latihan pemodelan atau dapatkah mereka berubah dari satu perhitungan ke yang berikutnya?
whuber
The p_is adalah tetap.
David B
Sehubungan dengan tanggapan saat ini, dapatkah Anda membagikan taksiran (a) jumlah p dan (b) jumlah kuadratnya? Nilai-nilai ini menentukan opsi Anda.
Whuber
@whuber: ini sangat bervariasi di setiap kasus. Ini bukan modul satu kali yang saya buat (sayangnya).
David B
@ David Tapi apakah Anda tidak bisa memberikan panduan, seperti rentang tipikal? Misalnya, jika jumlah p berkisar antara 1 dan 100, itu adalah informasi yang berguna dan menyarankan beberapa solusi yang efisien, tetapi jika bisa mencapai 10.000, itu bisa mengecualikan beberapa pendekatan.
Whuber

Jawaban:

24

Jika sering menyerupai Poisson , pernahkah Anda mencoba memperkirakannya dengan Poisson dengan parameter ?λ=pi

EDIT : Saya telah menemukan hasil teoritis untuk membenarkan ini, serta nama untuk distribusi : itu disebut distribusi binomial Poisson . Ketidaksetaraan Le Cam memberitahu Anda seberapa dekat distribusinya didekati dengan distribusi Poisson dengan parameter λ = Σ p i . Ini memberitahu Anda kualitas kira-kira ini diatur oleh jumlah kuadrat dari p i , untuk parafrase Steele (1994) . Jadi jika semua Anda p i s cukup kecil, seperti yang sekarang muncul mereka, itu harus baik pendekatan yang cukup.Yλ=pipipi

EDIT 2 : Seberapa kecil 'cukup kecil'? Yah, itu tergantung seberapa baik Anda membutuhkan perkiraan untuk menjadi! The artikel Wikipedia pada teorema Le Cam memberikan bentuk yang tepat dari hasil yang saya sebut di atas: jumlah dari perbedaan absolut antara fungsi massa probabilitas (PMF) dari dan PMF di atas Poisson distribusi tidak lebih dari dua kali jumlah tersebut dari kotak p i s. Hasil lain dari Le Cam (1960) mungkin lebih mudah untuk digunakan: jumlah ini juga tidak lebih dari 18 kali terbesar p i . Ada beberapa hasil seperti itu ... lihat Serfling (1978) untuk satu ulasan.Ypipi

onestop
sumber
1
+1 Bukan ide yang buruk. Kemungkinan campuran Poissons kecil akan melakukan pekerjaan dengan baik, tergantung pada bagaimana pertanyaannya diklarifikasi.
whuber
1
Saya memang berpikir untuk menyarankan distribusi binomial negatif, yang muncul sebagai campuran Gamma-Poisson, tetapi itu memiliki varian lebih besar dari rata-rata, sedangkan masalah ini memiliki varian lebih kecil dari rata-rata. Berdasarkan itu, saya tidak yakin apakah campuran Poissons akan bekerja, karena pasti campuran tersebut akan memiliki varian lebih besar dari rata-rata ??
onestop
@onestop Di mana dikatakan bahwa varians kurang dari rata-rata? Saya melewatkan pernyataan itu.
Whuber
Maaf ya, itu agak samar tapi komentar ini tidak memungkinkan banyak penjelasan. mpiktas ini adalah varians, yang kurang dari rata-rata, Σ p i . Hanya sedikit kurang jika p i 's yang rata-rata meskipun sangat kecil, sehingga standar Poisson mungkin baik cukup approx. Mungkin saya harus memperluas jawaban saya di atas .. tetapi kemudian utas percakapan menjadi membingungkan. Bn=pi(1pi)pipi
onestop
Apa yang Anda maksud dengan ? Bagaimana cara mendapatkan X i nilai-nilai? XiXi
David B
11

Saya menemukan pertanyaan Anda saat mencari solusi untuk masalah ini. Saya tidak puas dengan jawaban di sini, tetapi saya pikir ada solusi yang cukup sederhana yang memberi Anda distribusi yang tepat, dan sangat mudah ditelusuri.

Distribusi jumlah dua variabel acak diskrit adalah konvolusi kepadatannya. Jadi jika Anda memiliki mana Anda tahu P ( X ) dan P ( Y ) maka Anda dapat menghitung:Z=X+YP(X)P(Y)

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

(Tentu saja untuk variabel acak Bernoulli Anda tidak perlu pergi cukup hingga tak terbatas.)

Anda dapat menggunakan ini untuk menemukan distribusi tepat jumlah RV Anda. Pertama jumlah dua RV bersama-sama dengan menggabungkan PDF mereka (misalnya [0,3, 0,7] * [0,6, 0,4] = [0,18, 0,54, 0,28]). Kemudian gabungkan distribusi baru itu dengan Bernoulli PDF Anda berikutnya (mis. [0,18, 0,54, 0,28] * [0,5, 0,5] = [0,09, 0,36, 0,41, 0,14]). Terus ulangi ini sampai semua RV telah ditambahkan. Dan voila, vektor yang dihasilkan adalah PDF yang tepat dari jumlah semua variabel Anda.

Saya telah memverifikasi dengan simulasi bahwa ini menghasilkan hasil yang benar. Itu tidak bergantung pada asumsi asimptotik, dan tidak memiliki persyaratan bahwa prob Bernoulli kecil.

Mungkin juga ada beberapa cara untuk melakukan ini lebih efisien daripada lilitan yang berulang, tetapi saya belum memikirkannya secara mendalam. Saya harap ini membantu seseorang!

alex
sumber
2
Sudahkah Anda mencoba ini dengan variabel 40K ?? (Saya ingin tahu berapa jam atau hari perhitungan yang dibutuhkan ...)
Whuber
5
(+1) Saya menemukan cara untuk membuat ide ini berfungsi. Ini membutuhkan dua teknik: pertama, gunakan FFT untuk konvolusi; kedua, jangan lakukan secara berurutan, tetapi bagilah dan taklukkan: lakukan secara berpasangan, kemudian lakukan hasilnya dalam pasangan terputus-putus, dll. Algoritme sekarang berskala sebagai daripada O ( n 2 ) untuk n probabilitas. Misalnya, Mathematica dapat menghitung seluruh distribusi untuk 40.000 probabilitas hanya dalam 0,4 detik. (1.000.000 dihitung dalam 10,5 detik.) Saya akan memberikan kode dalam komentar tindak lanjut. O(nlogn)O(n2)n
whuber
7
Ini kode 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 &]]] Untuk menerapkannya, lakukan sesuatu seperti p = RandomReal[{0, 1}, 40000]; pp = multinomial[p];. Ini menciptakan probabilitas pdan kemudian menghitung distribusi yang tepat pp. NB Ketika rerata ptidak ekstrim, distribusinya sangat dekat dengan normal: yang mengarah pada algoritma yang jauh lebih cepat.
whuber
9

@onestop memberikan referensi yang bagus. Artikel Wikipedia tentang distribusi binomial Poisson memberikan rumus rekursif untuk menghitung distribusi probabilitas yang tepat; membutuhkan upaya . Sayangnya, ini adalah jumlah yang bolak-balik, sehingga secara numerik tidak stabil: sia-sia untuk melakukan perhitungan ini dengan aritmatika floating point. Untungnya, ketika p i kecil, Anda hanya perlu menghitung sejumlah kecil probabilitas, sehingga upaya ini benar-benar sebanding dengan O ( n log ( μ i p i ) ) . Ketelitian yang diperlukan untuk melakukan perhitungan dengan aritmatika rasional (O(n2)piO(nlog(ipi))yaitu, tepat,sehingga ketidakstabilan angka tidak menjadi masalah) tumbuh cukup lambat sehingga waktu keseluruhan mungkin masih sekitar . Itu layak.O(n2)

Sebagai tes, saya membuat array probabilitas untuk berbagai nilai n hingga n = 2 16 , yang merupakan ukuran dari masalah ini. Untuk nilai-nilai kecil n (hingga n = 2 12 ) waktu untuk perhitungan probabilitas yang tepat adalah dalam detik dan diskalakan secara kuadrat, jadi saya memberanikan perhitungan untuk n = 2 16pi=1/(i+1)nn=216nn=212n=216hingga tiga SD di atas rata-rata (probabilitas untuk 0, 1, ..., 22 keberhasilan). Butuh 80 menit (dengan Mathematica 8), sesuai dengan perkiraan waktu. (Probabilitas yang dihasilkan adalah pecahan yang pembilang dan penyebutnya masing-masing memiliki sekitar 75.000 digit!) Ini menunjukkan perhitungan dapat dilakukan.

Alternatifnya adalah menjalankan simulasi panjang (sejuta percobaan harus dilakukan). Itu hanya harus dilakukan sekali, karena tidak berubah.pi

whuber
sumber
9

(Karena pendekatan ini tidak tergantung pada solusi lain yang dipasang, termasuk yang telah saya posting, saya menawarkannya sebagai tanggapan terpisah).

Anda dapat menghitung distribusi tepat dalam hitungan detik (atau kurang) asalkan jumlah pnya kecil.

Kami telah melihat saran bahwa distribusi mungkin kira-kira Gaussian (dalam beberapa skenario) atau Poisson (dalam skenario lain). Either way, kita tahu rerata adalah jumlah dari p i dan variansnya σ 2 adalah jumlah dari p i ( 1 - p i ) . Oleh karena itu distribusi akan terkonsentrasi dalam beberapa standar deviasi dari rata-rata, katakanlah z SD dengan z antara 4 dan 6 atau sekitar itu. Karena itu kita hanya perlu menghitung probabilitas bahwa jumlah X sama dengan (bilangan bulat) k untuk k = μμpiσ2pi(1pi)zzXk hingga k = μ + z σ . Ketika sebagian besar p i kecil, σ 2 kira-kira sama dengan (tetapi sedikit kurang dari) μ , jadi untuk menjadi konservatif kita dapat melakukan perhitungan untuk k dalam interval [ μ - z k=μzσk=μ+zσpiσ2μk[μzμ,μ+zμ]. For example, when the sum of the pi equals 9 and choosing z=6 in order to cover the tails well, we would need the computation to cover k in [969,9+69] = [0,27], which is just 28 values.

The distribution is computed recursively. Let fi be the distribution of the sum of the first i of these Bernoulli variables. For any j from 0 through i+1, the sum of the first i+1 variables can equal j in two mutually exclusive ways: the sum of the first i variables equals j and the i+1st is 0 or else the sum of the first i variables equals j1 and the i+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
sumber
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
sumber
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<.
mpiktas
@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.


sumber
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
sumber
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
sumber
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
sumber
1
+1 Thank you for the useful quantitative information about the approximation bounds. Welcome to our site!
whuber