# Which probability distribution function applies to my problem?

OK the problem came from an AI course that I attended this year. You can find the problem description here

In short we have a number of birds in the sky that we need to shoot down.

In the main problem there are 90 white birds, 23 yellow, 22 blue, 22 green, 22 red and 1 black for a total of 180 birds.

My question is: How can I find the probability that I will have killed at least one of each species after X kills?

I know it is somehow connected to the hypergeometric distribution because that can tell you what is, for example, the probability that you will have hit the black bird after X number of kills, since we can treat each kill as sampling without replacement. But how do I do this for all the species?

I'm asking this because one of my strategies required at least one killed bird from each species to use as example data, and I'm not sure how I could get a number that would say:

"After X kills you have P probability that you have hit at least one bird from each species".

• The multivariate hypergometric distribution ( en.wikipedia.org/wiki/…) is what you're looking for – MånsT Oct 24 '12 at 11:31

If you want to find the probability that a uniformly random subset of size $X$ contains at least one of each species, then instead of using the multivariate hypergeometric distribution, I think you want to use inclusion-exclusion.

To be concrete, assume that $X = 50.$ The multivariate hypergeometric distribution tells you the probability that in a random subset of size $50$, you get exactly $12$ white, $2$ yellow, $5$ blue, $11$ green, $19$ red, and $1$ black bird. The problem is that you need to add these up over all possibilities with at least one of each type, and there are $137,083$ terms to add up. You could add up the complementary terms and subtract from $1$, but there are $244,692$ of those.

It is much simpler to use inclusion-exclusion. For each subset $S$ of the species, it is easy to calculate the probability $p(S)$ that no birds of that species were included in the $50$ birds. If $c(S)$ is the number of birds of species in $S$, then $p(S) = {180-c(S) \choose 50}/{180\choose 50}$. By inclusion-exclusion, the probability that you get at least one of each is

$$\sum_S (-1)^{|S|} p(S) = \sum_S (-1)^{|S|} \frac{180 - c(S) \choose X}{180 \choose X}.$$

There are $6$ species, so this sum is over the $2^6 = 64$ subsets of the set of species. That's a lot smaller than $137,083$, and the terms are simpler.

• +1 Very good! It looks to me that we had the same idea, you finished your answer while I was still working on mine. – gui11aume Oct 24 '12 at 14:53
• Well I will have to study this a bit more to fully understand it, but after reading the wiki pages on both subjects I think I have a better understanding. Thanks! – Bar Oct 25 '12 at 18:46

The annoying part here is that you have to sum the probabilities of the multivariate hypergeometric over all the 6-tuples that add up to $X$ and for which all 6 numbers are positive.

It is tempting to use an approximation which is much easier to calculate and that relies on the Poincaré formula

$$P(A_1 \cup ... \cup A_6) = P(A_1) + ... + P(A_6) - P(A_1 \cap A_2) - ...$$

in which all subsequent terms combine 2, 3, 4, 5 and 6 events. If we neglect those terms we end up with

$$P(A_1 \cup ... \cup A_6) \approx P(A_1) + ... + P(A_6)$$

The event $A_i$ "bird of color $i$ was not drawn" has probability $(1-\pi_i)^X$, where $\pi$ is the proportion of bird of color $i$. In the end you can approximate your probability by

$$(1-90/180)^X + (1-23/180)^X + 3 \times (1-22/180)^X + (1-1/180)^X.$$

However, this approximation is bad because it is higher than 1 for small $X$ and very far from 1 for $X=180$. As shown in this numeric example in R

X <- seq(from=10, to=180, by=10)
(1-90/180)^X + (1-23/180)^X + 3 * (1-22/180)^X + (1-1/180)^X


I think you are left with Monte Carlo simulations or with painful term-by-term computations.

Here is quick way to get it by simulation (takes about 1 second on my laptop).

birds <- c(rep(1,90), rep(2, 23), rep(3:5, each=22), 6)
X <- 20
got_them_all <- rep(NA, 10000)
set.seed(123)
for (i in 1:10000) {
got_them_all[i] <- length(unique(sample(x=birds, size=X, replace=FALSE))) == 6
}
mean(got_them_all) # answer: 0.0773


And here is a way that uses the full Poincaré formula

pi <- c(90, 23, rep(22,3), 1)/180
X <- 20
p <- 0
for (i in 1:6) {
p <- p + (-1)^(i+1) * sum((1-colSums(combn(pi,i)))^X)
}
1 - p # answer: 0.07403239