# Probability of drawing all balls in a non-uniformly distributed population

I am wondering if there is a way of calculating the following:

I have a bag with 5 balls numbered from 1 through 5 that are going to be drawn. Due to different weights, the balls have different probabilities of being drawn:

$P(Ball1) = 0.4$
$P(Ball2) = 0.3$
$P(Ball3) = 0.1$
$P(Ball4) = 0.1$
$P(Ball5) = 0.1$

There will be three draws. After each draw the balls are put back into the bag:

• Draw 1: Three balls
• Draw 2: Three balls
• Draw 3: Two balls

Now, is it possible to calculate the probability that each ball will have been chosen at least once along the three draws?

Disclaimer: I am not sure if the title is totally correct, but I think that my problem is a variant of Expected number of uniques in a non-uniformly distributed population - and therefore the similar title. Correct the title as you think it may be more correct.

My question is how to calculate this: $Pr(Ball1 > 0, Ball2 >0, Ball3 > 0, Ball4 >0, Ball5 > 0)$ given 3 + 3 + 2 balls drawn.

Edit: The balls that are drawn (3 + 3 + 2) are always drawn at once which would imply that the Fisher's noncentral hypergeometric distribution is what I am looking for. Note that if the drawing would happen ball by ball, Wallenius' noncentral hypergeometric distribution, would be the distribution of choice.

• I think the Multinomial distribution might help you...see here: en.wikipedia.org/wiki/Multinomial_distribution – Sam Livingstone Nov 15 '13 at 15:26
• @SamLivingstone, thank you. I have added the corresponding tag and am reading up. My problem seems similar to the Example in the Wikipedia entry. – dmeu Nov 15 '13 at 15:41
• Presumably the balls are only replaced after each set of 3. Otherwise you might just have drawn 8 with replacement. If so, you are going to have to decide whether the interation between drawing 3 at a time and the weights uses Wallenius' or Fisher's drawing method. – Henry Nov 15 '13 at 16:24
• Thank you @Henry, I was unaware of the two alternatives. I'd choose the Fisher's drawing method (suppose the 3 balls are being grabbed at once). – dmeu Nov 15 '13 at 16:29
• @dmeu: Fine (that means the drawn balls get put back unobserved if two or four come out and the draw repeated). There are $1000$ possible draws of 3 then 3 then 2, so you just have to work out the probability of each, and see which cover all five balls – Henry Nov 15 '13 at 16:43

There is no nice formula but there is a nice algorithm. The multiplication of polynomial generating functions mirrors the sequential drawing process. As such, computing with generating functions essentially performs an exhaustive enumeration of the possibilities. However, algebraic properties can be exploited to simplify the intermediate results (partial products), leading to a considerably streamlined calculation. These algebraic properties (formally: an ideal of polynomials) encode the two facts that (a) we are not concerned with how many of each ball is drawn, only which balls are drawn; and (b) it suffices to focus on the cases where not all balls are drawn.

This reply shows how these ideas can be used efficiently to compute an answer to the question--as well as to its obvious generalizations to more balls and different numbers and sizes of draws. It is followed by analyses, including a simulation, demonstrating the correctness of the answer. A final section discusses further issues that might be of interest, such as applying this approach to the simpler (unweighted) case.

To draw $k$ balls (with weighted probabilities) from an urn with $n$ balls we draw the first (but do not replace it), then the second (without replacement), and so on. Having observed the result, we replace all balls. Therefore the chance of seeing a sequence of balls with probabilities $p_1, p_2, \ldots, p_k$ equals

$$\Pr(p_1, p_2, \ldots, p_k) = p_1\frac{p_2}{1-p_1}\cdots\frac{p_k}{1-(p_1+p_2+\cdots+p_{k-1})}$$

and the chance of seeing a combination (unordered sequence) of balls with those probabilities is the sum over all permutations,

$$\Pr(\{p_1, p_2, \ldots, p_k\}) = \sum_{\sigma\in\mathcal{S}_k} \Pr(p_{\sigma(1)}, p_{\sigma(2)}, \ldots, p_{\sigma(k)}).$$

A simple and elegant way to handle these $\binom{n}{k}$ numbers is with a generating polynomial,

$$m_k(x_1, x_2, \ldots, x_n) = \sum_{i_1\lt i_2\lt \cdots \lt i_k}\Pr(\{p_{i_1}, p_{i_2}, \ldots, p_{i_k}\})x_{i_1}x_{i_2}\cdots x_{i_k}.$$

For instance, with probabilities $(4/10, 3/10, 1/10, 1/10, 1/10),$ we may compute

$$m_1(x_1, x_2, x_3, x_4, x_5) = \frac{4 x_1}{10}+\frac{3 x_2}{10}+\frac{x_3}{10}+\frac{x_4}{10}+\frac{x_5}{10},$$

$$m_2(x_1, x_2, x_3, x_4, x_5) = \frac{13 x_1 x_2}{35}+\frac{8 x_3 x_2}{105}+\frac{8 x_4 x_2}{105}+\cdots+\frac{x_3 x_5}{45}+\frac{x_4 x_5}{45},$$

and

$$m_3(x_1, x_2, x_3, x_4, x_5) =\frac{76}{315} x_1 x_2 x_3+\frac{76}{315} x_1 x_2 x_4\cdots+\frac{17}{504} x_2 x_4 x_5+\frac{1}{120} x_3 x_4 x_5.$$

Writing $x = (x_1, x_2, \ldots, x_n),$ the generating polynomial for the total outcome of a sequence of independent draws of sizes $k_1, k_2, \ldots, k_d$ is simply the product $m_{k_1}(x)m_{k_2}(x)\cdots m_{k_d}(x).$ Moreover, all we are interested in is whether all the $x_i$ appear in the outcome, not how many times each does. This allows the computation of the product to be simplified: at each multiplication we may throw away any monomial term that is a multiple of $x_1x_2\cdots x_n$ and replace any power of any $x_i$ by $x_i$ itself (which can be done by equating each $x_i^2$ with $x_i$ itself). The result will encode the chances of not obtaining all $n$ balls among these $d$ draws: the coefficient of each monomial $x_{i_1}x_{i_2}\cdots x_{i_s}$ will be the chance of obtaining exactly the balls $\{i_1, i_2, \ldots, i_s\}$. Setting all the $x_i=1$ will give the chance that not all balls are drawn.

The calculations are easiest to perform with a computer algebra system. In Mathematica a solution looks like this:

1. Compute the coefficients of the generating polynomials recursively:

prob[{j_, k___}, p_List, x_] := Block[{n = Total@p, q = p},
q[[j]] = 0;
prob[{k}, q, x] p[[j]] Subscript[x, j] / n] ;
prob[{}, _, _] := 1

2. Sum over all permutations to obtain the generating polynomials:

m[k_Integer, p_, x_] := m[k, p, x] =
Sum[prob[y, p, x], {y, Permutations[Range[Length[p]], {k}]}]

3. Specify the ideal generated by $x_1x_2\cdots x_n$ and the $x_i^2-x_i$:

ideal[p_, x_] :=
With[{n = Length[p]}, {Product[Subscript[x, i], {i, 1, n}],
Subscript[x, #]^2 - Subscript[x, #] & /@ Range[n]} // Flatten]

4. Multiply the generating polynomials, reducing modulo the ideal at each step:

Block[{weights = {4, 3, 1, 1, 1}, draws = {3, 3, 2}, x, i, gf},
i = ideal[weights, x];
gf = Fold[PolynomialMod[#1 m[#2, weights, x], i] &, 1, draws];

5. Set all $x_i$ to $1$:

 gf /. (Subscript[x, #] -> 1 & /@ Range[Length[weights]])]


The answer takes $0.02$ seconds to compute and equals $\frac{74344589}{111132000}$. This means the chance of obtaining all five balls is $1 -\frac{74344589}{111132000} = \frac{36787411}{111132000} \approx 0.33102447.$

### Discussion

Simulation

Whenever possible, it is a good idea to double-check complicated or tricky probability calculations with a simulation. Here is a quick simulation (for arbitrary weights and arbitrary numbers of draws) in R, as applied to the original problem:

weights <- c(4,3,1,1,1)
draws <- c(3,3,2)
n.iter <- 10^4
set.seed(17)
p <- mean(replicate(n.iter, {
length(table(unlist(sapply(draws,
function(d) sample.int(5, d, prob=weights))))) >= length(weights)
}))
cat("Estimate =", p,
"( Z =", (p - 36787411/111132000) / sqrt(p*(1-p)/n.iter), ")")


Estimate = 0.335 ( Z = 0.8422911 )

The output reports the proportion of iterations in which all the balls were observed, followed by a Z-statistic relative to the previously-computed answer. Because the Z-statistic is small, this is evidence that the computed answer is correct (up to a few times the standard error, which in this example is $0.005$).

Application to simpler problems

One case in which there is a formula (of sorts) is when there are no weights. The Mathematica code, applied to a weights vector of {1,1,1,1,1}, outputs $\frac{9}{20}$, yielding an answer of $1 - \frac{9}{20} = \frac{11}{20} = 0.55.$ Is it right?

The symmetry of this problem--all balls have equal probabilities--makes an independent calculation of the answer feasible. After drawing the first three balls, name them $x_1, x_2,$ and $x_3$, and let the remaining balls be $x_4$ and $x_5$. Replace all balls and draw another three. There are two relevant events to track: whether $3$ or $4$ balls have been observed. The chance of $3$ is the chance that the same $3$ balls were drawn a second time; this is $1 / \binom{5}{3} = 1/10.$ The chance of observing $4$ distinct balls is that two of the new balls were from $\{x_1, x_2, x_3\}$ and the third from $\{x_4, x_5\}$: this equals $\frac{\binom{3}{2}\binom{2}{1}}{\binom{5}{3}} = \frac{3\times 2}{10} = \frac{6}{10}.$

In each of these cases we compute the chance that the final draw of $2$ balls will not complete the set of five. If only $3$ balls have been observed, that chance is $9/10$. If $4$ balls have been observed, the chance is $\frac{\binom{4}{2}}{\binom{5}{2}} = \frac{6}{10}.$ Whence the chance of not completing the set is

$$\frac{1}{10} \frac{9}{10} + \frac{6}{10}\frac{6}{10} = \frac{9+36}{100} = \frac{45}{100} = \frac{9}{20},$$

exactly as reported by the Mathematica code.

Closer examination of the polynomials

How the code works is easy to see in this case, because the polynomials $m_k$ are multiples of the elementary symmetric functions:

$$10 m_2(x) = x_1 x_2+x_3 x_2+x_4 x_2+x_5 x_2+x_1 x_3+x_1 x_4+x_3 x_4+x_1 x_5+x_3 x_5+x_4 x_5$$

$$10 m_3(x) = x_1 x_2 x_3+x_1 x_4 x_3+x_2 x_4 x_3+\cdots+x_1 x_2 x_5+x_1 x_4 x_5+x_2 x_4 x_5.$$

These are constructed by taking a single product, such as $x_1x_2$, and adding to it all monomials attainable by permuting their subscripts (among the set $\{1,2,3,4,5\}$). In effect, each monomial has only one kind of term plus its permuted versions.

Let's compute their products modulo the ideal $\mathcal{I}$ generated by $\{x_1x_2x_3x_4x_5, x_1^2-x_1, x_2^2-x_2, x_3^2-x_3, x_4^2-x_4, x_5^2-x_5\}.$ To do so we only have to look at a limited number of products and then apply all permutations. For instance, in computing $(10 m_3(x))^2$ we find one term like $(x_1x_2x_3)^2 \equiv x_1x_2x_3$ mod $\mathcal{I}$--together with all its permutations--and terms that are congruent to $x_1x_2x_3x_4$ mod $\mathcal{I}$--together with all their permutations. It's straightforward to count that there are $12$ each of the latter, whence

$$(10 m_3(x))^2 \equiv\sum_{\sigma\in \mathcal{S}_5}\left(x_{\sigma(1)}x_{\sigma(2)}x_{\sigma(3)} + 12x_{\sigma(1)}x_{\sigma(2)}x_{\sigma(3)}x_{\sigma(4)}\right)\quad\text{mod }\mathcal{I}.$$

Multiplication by $10m_2(x)$ modulo $\mathcal{I}$ is similarly easy. The result is

$$(10 m_3(x))^2(10 m_2(x)) \equiv \sum_{\sigma\in \mathcal{S}_5}\left(3x_{\sigma(1)}x_{\sigma(2)}x_{\sigma(3)} + 84x_{\sigma(1)}x_{\sigma(2)}x_{\sigma(3)}x_{\sigma(4)}\right)\quad\text{mod }\mathcal{I}.$$

This says there are $\binom{5}{3}\times 3$ chances in $1000$ of observing exactly three of the balls after independent draws of $3$, $3$, and $2$ balls and $\binom{5}{4}\times 84$ chances in $1000$ of observing exactly four of the balls. Their sum is $10\times 3 + 5\times 84 = 450$ chances out of $1000$.

The general method used in the weighted-probability case performs a similar set of computations.

Alternative assumptions

For alternative drawing procedures (see the comment thread for the question), merely change the definition of prob, which describes the chances of drawing each sequence of balls. The rest of the algorithm remains unchanged.

• Wow, thanks @whuber. I have worked it out also, but I guess not nearly as elegant as you. My calculation gave me p = $0.3993262$. I was using Fisher's noncentral hypergeometrical distribution (because I am grabbing the balls in a draw at once). You did use Wallenius, right? I will simplify the code and put it as another answer so you can review it. It is in R (and I do not understand much Mathematica). – dmeu Nov 18 '13 at 10:37
• I have tried it with the Wallanius distribution, and obtained the same probability as @whuber, but in my specific example would be the wrong one. May I ask how your algorithm behaves with bigger numbers? Let's say that we want to draw all 100 different balls by drawing 60, 35 and 20 balls in three draws? For my algorithm these numbers are already nonviable. – dmeu Nov 18 '13 at 12:44
• That completely changes the question, because it is unrealistic to expect to compute an exact answer in the lifetime of the universe. (Just computing the probability distribution for $60$ out of $100$ balls with all unequal weights requires over $10^{28}$ multiplications; the full probability distribution for $60, 35$ and $20$ balls needs almost $10^{76}$ multiplications.) You are effectively asking for some sort of approximation, which will require different techniques (as hinted at in the Wikipedia articles referenced by @Henry in a comment to your question). – whuber Nov 18 '13 at 16:08
• Ok, thanks for your answer. I am learning step by step. I will look into it to see what I can figure out and if necessary open a new thread referencing this one. – dmeu Nov 18 '13 at 16:17
• I was troubled with something similar recently, but was able to run it to sample sizes up to a few thousands using Wallenius' Markov approximation. More info here and in his original PhD thesis and Agner Fog's work. – ciri Jan 30 '14 at 18:11

Thus, with lots of help in the comments of the question, I have worked out the following solution in R. I think the result is correct, but am glad if anybody could review it.

The algorithm from @whuber may be less costly to calculate. As soon as the amount of different balls and N of drawn balls rise, the the combn() function will fail due to too many combinations.

The result of the example in the question is p = 0.3993262 with Fisher's noncentral hypergeometric destribution and p = 0.3310244 with Wallenius' noncentral hypergeometric distribution.

Define the variables

<!-- language: lang-R -->
weights = c(0.4, 0.3, 0.1, 0.1, 0.1)
names = c("draw1","draw2","draw3")
drawN = c(3, 3, 2)


Create empty list to store data:

 <!-- language: lang-R -->

N = length(weights)
draws = list()
probs = list()
groups = list()


Calculate the probabilities for each possible draw of the three draws: 3 times 10 possibilities. As pointed out by @henry in the comments, I am using the Fisher's non-central hypergeometric distribution: dMFNCHypergeo()

  <!-- language: lang-r -->
for (i in seq(1,length(names))) {
name <- names[i]
d <- data.frame(combn(seq(1,N),drawN[i]))
colnames(d) <- paste(name,colnames(d),sep=".")
draws <- c(draws,d)
groups[[name]] <- colnames(d)

urn <- rep(1,N)
for (col in colnames(d)) {
thisdraw <- rep(0,N)
thisdraw[d[,col]] <- 1
probs[[col]] <- dMFNCHypergeo(thisdraw, urn, sum(thisdraw), weights)
}
}


Combine the different draws and sum the wanted probabilities: 10 * 10 * 10 = 1000

  <!-- language: lang-R -->
groupcombs <- as.matrix(expand.grid(groups))
groupprobs <- unlist(probs)

psum = 0;
is <- c()
for (i in seq(1,dim(groupcombs)[1])) {
d <- draws[groupcombs[i,]]
if (length(unique(unlist(d))) == N) {
p <- prod(groupprobs[names(d)])
psum <- psum + p
}
}
print(psum)

• Why would Fisher's distribution be appropriate? It does not seem to correspond to your description of the ball-drawing process. Perhaps you could edit the question to clarify which process is relevant. – whuber Nov 18 '13 at 15:30
• Yes, that Fisher's distribution is what I was looking for only became evident to me in the discussions in the comments of the question. I will put it in the question. – dmeu Nov 18 '13 at 16:15