Distribusi probabilitas untuk berbagai probabilitas

36

Jika saya ingin mendapatkan probabilitas 9 keberhasilan dalam 16 percobaan dengan setiap percobaan memiliki probabilitas 0,6 saya bisa menggunakan distribusi binomial. Apa yang bisa saya gunakan jika masing-masing dari 16 percobaan memiliki probabilitas keberhasilan yang berbeda?

Greg
sumber
1
@whuber Dalam penjelasan Anda tentang perkiraan Normal, perhitungan rata-rata dan standar deviasi berbeda dengan deskripsi di Wikipedia. Di Wiki, rata-rata adalah np, dan standar deviasi adalah np (1-p). Jadi, dalam masalah ini, untuk perkiraan Normal dari berbagai kemungkinan keberhasilan dalam distribusi binomial, rerata adalah p1 + p2 + p3 + p4 + p5 + ... + pi, dan variansnya adalah p1 (1-p1) + p2 ( 1-p2) + ... + pi (1-pi). Apakah saya benar?
David
1
Lihat Wikipedia di distribusi binomial Poisson . Juga istilah pencarian yang menghasilkan beberapa klik di sini.
Glen_b -Reinstate Monica
@David Bila semua pi sama dengan nilai umum p , maka p1+p2++pn=np dan p1(1p1)++pn(1pn)=np(1p) , menunjukkan bahwa deskripsi Wikipedia yang Anda rujuk hanyalah kasus khusus.
whuber

Jawaban:

22

Ini adalah jumlah dari 16 uji coba Binomial (mungkin independen). Asumsi kemerdekaan memungkinkan kita untuk melipatgandakan probabilitas. Dimana, setelah dua percobaan dengan probabilitas p1 dan p2 keberhasilan peluang keberhasilan pada kedua percobaan adalah p1p2 , peluang tidak ada keberhasilan adalah (1p1)(1p2) , dan peluang untuk satu kesuksesan adalah p1(1p2)+(1p1)p2 . Ungkapan terakhir itu berutang validitasnya pada kenyataan bahwa dua cara untuk mendapatkan satu kesuksesan itu sama-sama eksklusif: paling-paling satu di antaranya bisa benar-benar terjadi. Itu berarti probabilitas merekabertambah.

Dengan menggunakan dua aturan ini - probabilitas independen bertambah banyak dan yang eksklusif tambahkan - Anda dapat menentukan jawabannya, katakanlah, 16 uji coba dengan probabilitas p1,,p16 . Untuk melakukannya, Anda harus memperhitungkan semua cara untuk memperoleh setiap keberhasilan yang diberikan (seperti 9). Ada (169)=11440cara untuk mencapai 9 keberhasilan. Salah satunya, misalnya, terjadi ketika percobaan 1, 2, 4, 5, 6, 11, 12, 14, dan 15 berhasil dan yang lainnya gagal. Keberhasilan memiliki probabilitasp1,p2,p4,p5,p6,p11,p12,p14,danp15dan kegagalan memiliki probabilitas1p3,1p7,,1p13,1p16 . Mengalikan 16 angka ini memberi peluangurutan hasil tertentu. Menjumlahkan angka ini bersama dengan 11,439 angka yang tersisa memberikan jawabannya.

Tentu saja Anda akan menggunakan komputer.

Dengan lebih dari 16 percobaan, ada kebutuhan untuk memperkirakan distribusi. Asalkan tidak ada probabilitas pi dan 1pi menjadi terlalu kecil, perkiraan Normal cenderung berfungsi dengan baik. Dengan metode ini Anda mencatat bahwa ekspektasi jumlah penjumlahan n percobaan adalah μ=p1+p2++pn dan (karena uji coba independen) variansnya adalah σ2=p1(1p1)+p2(1p2)++pn(1pn) . Anda kemudian berpura-pura distribusi jumlah adalah Normal dengan meanμ dan simpangan bakuσ . Jawabannya cenderung baik untuk menghitung probabilitas yang sesuai dengan proporsi keberhasilan yang berbeda dariμ dengan tidak lebih dari beberapa kelipatanσ . Saatn tumbuh besar, perkiraan ini menjadi lebih akurat dan bekerja untuk kelipatanσ jauh lebih besar dariμ .

whuber
sumber
9
Ilmuwan komputer menyebutnya "uji coba Poisson" untuk membedakannya dari uji coba Bernoulli. Selain perkiraan Central Limit Theorem, ada juga batas ekor yang baik. Ini satu. Pencarian Google pada "batas Chernoff untuk uji coba Poisson" akan muncul hasil yang mungkin Anda temukan dalam perawatan CS khas.
kardinal
@ Cardinal Nomenklatur itu menarik. Ini akan berlaku untuk sangat kecil pi, tetapi sebaliknya keliru, karena distribusi sebaliknya tidak didekati dengan baik oleh distribusi Poisson. (Ada diskusi lain di CV tentang pertanyaan ini, di mana "16" diganti dengan 10.000 dan kami memang memeriksa probabilitas ekor, tetapi saya belum dapat menemukannya lagi.)
whuber
1
ya, saya setuju namanya. Saya merasa sedikit aneh ketika saya pertama kali menemukannya. Saya telah memberikannya di sini lebih sebagai istilah yang berguna untuk pencarian. Tampaknya para ilmuwan komputer sering mempertimbangkan probabilitas ini dalam berurusan dengan algoritma tertentu. Saya akan tertarik membaca pertanyaan lain itu jika Anda menemukannya. Apakah yang ini mungkin?
kardinal
2
@ cardinal benar bahwa kita "orang-orang CS" menyebutnya uji coba Poisson. sebenarnya untuk kasus ini, standar Chernoff-Hoeffding terikat akan memberi Anda persis batas yang diminta OP.
Suresh Venkatasubramanian
1
Sesuai komentar oleh @David kemarin, ada sesuatu yang salah dengan pernyataan rata-rata perkiraan Normal Anda sebagai Kami menjumlahkan 16 Bernoulli rvs, masing-masing dapat mengambil nilai 0 atau 1, sehingga jumlah akan memiliki domain dukungan dari 0 hingga 16, bukan antara 0 dan 1. Layak memeriksa sd Anda juga.
μ=(p1+p2++pn)/n
serigala
12

Salah satu alternatif untuk pendekatan normal @ whuber adalah dengan menggunakan probabilitas "mixing", atau model hierarkis. Ini akan berlaku ketika serupa dalam beberapa cara, dan Anda dapat memodelkan ini dengan distribusi probabilitas p iD i s t ( θ ) dengan fungsi kerapatan g ( p | θ ) yang diindeks oleh beberapa parameter θ . Anda mendapatkan persamaan integral:pipiDist(θ)g(p|θ)θ

Pr(s=9|n=16,θ)=(169)01p9(1p)7g(p|θ)dp

Probabilitas binomial berasal dari pengaturan , perkiraan normal berasal dari (saya pikir) pengaturan g ( pg(p|θ)=δ(pθ)(denganμdanσsebagaimana didefinisikan dalam jawaban @ whuber) dan kemudian mencatat "ekor" dari PDF ini jatuh tajam di sekitar puncak.g(p|θ)=g(p|μ,σ)=1σϕ(pμσ)μσ

Anda juga bisa menggunakan distribusi beta, yang akan mengarah ke bentuk analitik sederhana, dan yang tidak perlu menderita masalah "p kecil" yang dilakukan oleh perkiraan normal - karena beta cukup fleksibel. Menggunakan distribusi dengan α , β yang ditetapkan oleh solusi untuk persamaan berikut (ini adalah perkiraan "mimimum KL divergence"):beta(α,β)α,β

ψ(β)-ψ(α+β)=1

ψ(α)ψ(α+β)=1ni=1nlog[pi]
ψ(β)ψ(α+β)=1ni=1nlog[1pi]

Di mana Adalah fungsi digamma - terkait erat dengan deret harmonik.ψ(.)

We get the "beta-binomial" compound distribution:

(169)1B(α,β)01p9+α1(1p)7+β1dp=(169)B(α+9,β+7)B(α,β)

This distribution converges towards a normal distribution in the case that @whuber points out - but should give reasonable answers for small n and skewed pi - but not for multimodal pi, as beta distribution only has one peak. But you can easily fix this, by simply using M beta distributions for the M modes. You break up the integral from 0<p<1 into M pieces so that each piece has a unique mode (and enough data to estimate parameters), and fit a beta distribution within each piece. then add up the results, noting that making the change of variables p=xLUL for L<x<U the beta integral transforms to:

B(α,β)=LU(xL)α1(Ux)β1(UL)α+β1dx
probabilityislogic
sumber
+1 This answer contains some interesting and clever suggestions. The last one looks particularly flexible and powerful.
whuber
Just to take something very simple and concrete, suppose (i) pi=i17 and (ii) pi=i/17, for i=1 to 16. What would be the solution to your α and β estimates, and thus your estimates for P(X=9) given n=16, as per the OP's problem?
wolfies
Great answer and proposal, especially the beta! It'd be cool to see this answer written in its general form with n and s.
pglpm
8

XiBernoulli(pi)

pgf=E[tXi]=1pi(1t)

S=i=1nXi denote the sum of n such independent random variables. Then, the pgf for the sum S of n=16 such variables is:

pgfS=E[tS]=E[tX1]E[tX2]E[tX16] (... by independence)=i=116(1pi(1t))

We seek P(S=9), which is:

19!d9pgfSdt9|t=0

ALL DONE. This produces the exact symbolic solution as a function of the pi. The answer is rather long to print on screen, but it is entirely tractable, and takes less than 1100th of a second to evaluate using Mathematica on my computer.

Examples

If pi=i17,i=1 to 16, then: P(S=9)=964794185433480818448661191875666868481=0.198268

If pi=i17,i=1 to 16, then: P(S=9)=0.000228613

More than 16 trials?

With more than 16 trials, there is no need to approximate the distribution. The above exact method works just as easily for examples with say n=50 or n=100. For instance, when n=50, it takes less than 110th of second to evaluate the entire pmf (i.e. at every value s=0,1,,50) using the code below.

Mathematica code

Given a vector of pi values, say:

n = 16;   pvals = Table[Subscript[p, i] -> i/(n+1), {i, n}];

... here is some Mathematica code to do everything required:

pgfS = Expand[ Product[1-(1-t)Subscript[p,i], {i, n}] /. pvals];
D[pgfS, {t, 9}]/9! /. t -> 0  // N

0.198268

To derive the entire pmf:

Table[D[pgfS, {t,s}]/s! /. t -> 0 // N, {s, 0, n}]

... or use the even neater and faster (thanks to a suggestion from Ray Koopman below):

CoefficientList[pgfS, t] // N

For an example with n=1000, it takes just 1 second to calculate pgfS, and then 0.002 seconds to derive the entire pmf using CoefficientList, so it is extremely efficient.

wolfies
sumber
1
It can be even simpler. With[{p = Range@16/17}, N@Coefficient[Times@@(1-p+p*t),t,9]] gives the probability of 9 successes, and With[{p = Range@16/17}, N@CoefficientList[Times@@(1-p+p*t),t]] gives the probabilities of 0,...,16 successes.
Ray Koopman
@RayKoopman That is cool. The Table for the p-values is intentional to allow for more general forms not suitable with Range. Your use of CoefficientList is very nice! I've added an Expand to the code above which speeds the direct approach up enormously. Even so, CoefficientList is even faster than a ParallelTable. It does not make much difference for n under 50 (both approaches take just a tiny fraction of a second either way to generate the entire pmf), but your CoefficientList will also be a real practical advantage when n is really large.
wolfies
5

@wolfies comment, and my attempt at a response to it revealed an important problem with my other answer, which I will discuss later.

Specific Case (n=16)

There is a fairly efficient way to code up the full distribution by using the "trick" of using base 2 (binary) numbers in the calculation. It only requires 4 lines of R code to get the full distribution of Y=i=1nZi where Pr(Zi=1)=pi. Basically, there are a total of 2n choices of the vector z=(z1,,zn) that the binary variables Zi could take. Now suppose we number each distinct choice from 1 up to 2n. This on its own is nothing special, but now suppose that we represent the "choice number" using base 2 arithmetic. Now take n=3 so I can write down all the choices so there are 23=8 choices. Then 1,2,3,4,5,6,7,8 in "ordinary numbers" becomes 1,10,11,100,101,110,111,1000 in "binary numbers". Now suppose we write these as four digit numbers, then we have 0001,0010,0011,0100,0101,0110,0111,1000. Now look at the last 3 digits of each number - 001 can be thought of as (Z1=0,Z2=0,Z3=1)Y=1, etc. Counting in binary form provides an efficient way to organise the summation. Fortunately, there is an R function which can do this binary conversion for us, called intToBits(x) and we convert the raw binary form into a numeric via as.numeric(intToBits(x)), then we will get a vector with 32 elements, each element being the digit of the base 2 version of our number (read from right to left, not left to right). Using this trick combined with some other R vectorisations, we can calculate the probability that y=9 in 4 lines of R code:

exact_calc <- function(y,p){
    n       <- length(p)
    z       <- t(matrix(as.numeric(intToBits(1:2^n)),ncol=2^n))[,1:n] #don't need columns n+1,...,32 as these are always 0
    pz      <- z%*%log(p/(1-p))+sum(log(1-p))
    ydist   <- rowsum(exp(pz),rowSums(z))
    return(ydist[y+1])
}

Plugging in the uniform case pi(1)=i17 and the sqrt root case pi(2)=i17 gives a full distribution for y as:

yPr(Y=y|pi=i17)Pr(Y=y|pi=i17)00.00000.055810.00000.178420.00030.265230.00260.243040.01390.153650.04910.071060.11810.024870.19830.006780.23530.001490.19830.0002100.11810.0000110.04910.0000120.01390.0000130.00260.0000140.00030.0000150.00000.0000160.00000.0000

So for the specific problem of y successes in 16 trials, the exact calculations are straight-forward. This also works for a number of probabilities up to about n=20 - beyond that you are likely to start to run into memory problems, and different computing tricks are needed.

Note that by applying my suggested "beta distribution" we get parameter estimates of α=β=1.3206 and this gives a probability estimate that is nearly uniform in y, giving an approximate value of pr(y=9)=0.06799117. This seems strange given that a density of a beta distribution with α=β=1.3206 closely approximates the histogram of the pi values. What went wrong?

General Case

I will now discuss the more general case, and why my simple beta approximation failed. Basically, by writing (y|n,p)Binom(n,p) and then mixing over p with another distribution pf(θ) is actually making an important assumption - that we can approximate the actual probability with a single binomial probability - the only problem that remains is which value of p to use. One way to see this is to use the mixing density which is discrete uniform over the actual pi. So we replace the beta distribution pBeta(a,b) with a discrete density of pi=116wiδ(ppi). Then using the mixing approximation can be expressed in words as choose a pi value with probability wi, and assume all bernoulli trials have this probability. Clearly, for such an approximation to work well, most of the pi values should be similar to each other. This basically means that for @wolfies uniform distribution of values, pi=i17 results in a woefully bad approximation when using the beta mixing distribution. This also explains why the approximation is much better for pi=i17 - they are less spread out.

The mixing then uses the observed pi to average over all possible choices of a single p. Now because "mixing" is like a weighted average, it cannot possibly do any better than using the single best p. So if the pi are sufficiently spread out, there can be no single p that could provide a good approximation to all pi.

One thing I did say in my other answer was that it may be better to use a mixture of beta distributions over a restricted range - but this still won't help here because this is still mixing over a single p. What makes more sense is split the interval (0,1) up into pieces and have a binomial within each piece. For example, we could choose (0,0.1,0.2,,0.9,1) as our splits and fit nine binomials within each 0.1 range of probability. Basically, within each split, we would fit a simple approximation, such as using a binomial with probability equal to the average of the pi in that range. If we make the intervals small enough, the approximation becomes arbitrarily good. But note that all this does is leave us with having to deal with a sum of indpendent binomial trials with different probabilities, instead of Bernoulli trials. However, the previous part to this answer showed that we can do the exact calculations provided that the number of binomials is sufficiently small, say 10-15 or so.

To extend the bernoulli-based answer to a binomial-based one, we simply "re-interpret" what the Zi variables are. We simply state that Zi=I(Xi>0) - this reduces to the original bernoulli-based Zi but now says which binomials the successes are coming from. So the case (Z1=0,Z2=0,Z3=1) now means that all the "successes" come from the third binomial, and none from the first two.

Note that this is still "exponential" in that the number of calculations is something like kg where g is the number of binomials, and k is the group size - so you have Yj=1gXj where XjBin(k,pj). But this is better than the 2gk that you'd be dealing with by using bernoulli random variables. For example, suppose we split the n=16 probabilities into g=4 groups with k=4 probabilities in each group. This gives 44=256 calculations, compared to 216=65536

By choosing g=10 groups, and noting that the limit was about n=20 which is about 107 cells, we can effectively use this method to increase the maximum n to n=50.

If we make a cruder approximation, by lowering g, we will increase the "feasible" size for n. g=5 means that you can have an effective n of about 125. Beyond this the normal approximation should be extremely accurate.

probabilityislogic
sumber
@momo - I think this is ok, as my answers are two different ways to approach the problem. This answer is not an edited version of my first one - it is just a different answer
probabilityislogic
1
For a solution in R that is extremely efficient and handles much, much larger values of n, please see stats.stackexchange.com/a/41263. For instance, it solved this problem for n=104, giving the full distribution, in under three seconds. (A comparable Mathematica 9 solution--see @wolfies' answer--also performs well for smaller n but could not complete the execution with such a large value of n.)
whuber
5

The (in general intractable) pmf is

Pr(S=k)=A{1,,n}|A|=k(iApi)(j{1,,n}A(1pj)).
R code:
p <- seq(1, 16) / 17
cat(p, "\n")
n <- length(p)
k <- 9
S <- seq(1, n)
A <- combn(S, k)
pr <- 0
for (i in 1:choose(n, k)) {
    pr <- pr + exp(sum(log(p[A[,i]])) + sum(log(1 - p[setdiff(S, A[,i])])))
}
cat("Pr(S = ", k, ") = ", pr, "\n", sep = "")

For the pi's used in wolfies answer, we have:

Pr(S = 9) = 0.1982677

When n grows, use a convolution.

Zen
sumber
1
Doing that with R code was really helpful. Some of us are more concrete-thinkers and it greatly helps to have an operational version of the generating function.
DWin
@DWin I provide efficient R code in the solution to the same problem (with different values of the pi) at stats.stackexchange.com/a/41263. The problem here is solved in 0.00012 seconds total computation time (estimated by solving it 1000 times) compared to 0.53 seconds (estimated by solving it once) for this R code and 0.00058 seconds using Wolfies' Mathematica code (estimated by solving it 1000 times).
whuber
So P(S=k) would follow a Poisson-Binomial Distribution.
fccoelho
+1 Very useful post in my attempt at answering this question. I was wondering if using logs is more of a cool mathematical formulation than a real need. I am not too concerned about running times...
Antoni Parellada