Probability distribution over items belonging to a subset of boxes

I am trying to find a probability distribution that would help me answer the following question:

I have N different boxes and K different items. For each item, I have a probability distribution over the N boxes. The probability distributions for each item are different and independent. Given a subset of M boxes (M <= N), what is the probability that all K items a) belong only to the subset of M boxes and b) each box in the subset has at least one item.

Thank you.

• Note that N < 100 and K > 1000. Approximation might be necessary in this case. – ondrejba Sep 10 '18 at 20:16
• There's little problem with large $K,$ but with $N \gt 30$ or so there are just too many calculations to perform. Unless there is some regular pattern to the probabilities, why not just simulate the answer? – whuber Sep 10 '18 at 20:22
• How can I simulate the answer? – ondrejba Sep 10 '18 at 20:43
• You fill the boxes according to the distributions given and note whether any are empty. Repeat a few thousand times. You can focus on the given subset, because your event assumes all the items are in it. – whuber Sep 10 '18 at 21:15
• I tested the simulation but it's too slow for my application. How would I go about calculating it if N < 30? Thanks! – ondrejba Sep 11 '18 at 17:51

Suppose you have a set of $k$ discrete probabilities $p_1,p_2,\ldots,p_k$ defined on the sample space $\mathcal S.$ Extend the notation to all subsets $\mathcal{A}\subset \mathcal{S}$ by defining

$$p_i(\mathcal{A}) = \sum_{\omega\in\mathcal{A}} p_i(\omega)$$

to be the probabilities of $\mathcal A.$ Let

$$p(\mathcal{A}) = \prod_{i=1}^k p_i(\mathcal{A})$$

be the probability of $\mathcal{A}$ simultaneously for all $k$ independent distributions.

For natural numbers $n=0, 1, 2, \ldots,$ write

$$\mathcal{S}_{(-n)} = \{\mathcal{A}\subset\mathcal{S}\mid |\mathcal S \setminus \mathcal A|=n\}$$

for the collection of subsets of $\mathcal S$ whose complements have $n$ elements.

The Inclusion-Exclusion Principle asserts that the chance of each element of $\mathcal S$ occurring is

$$\sum_{n=0}^{|\mathcal S|}(-1)^n \sum_{\mathcal{A}\in\mathcal{S}_{(-n)}}p(\mathcal A).\tag{*}$$

In general there is no algebraic simplification. The number of computations scales as $2^{|\mathcal {S}|},$ which (being worse than polynomial) is terrible. Don't try this with sets much larger than 30 or so! The alternation of addition and subtraction is even worse, due to its potential to lose precision in floating point arithmetic.

Oh, yes, one more thing. The original question concerns a subset $\mathcal A$ of the bins $\mathcal S.$ The chance all items are in $\mathcal A$ is the product of the $p_i(\mathcal A),$ because their dispositions are independent. Replace the probabilities by their truncated versions, obtained by dividing $p_i(\omega)$ by $p_i(\mathcal A)$ for each $\omega\in\mathcal A,$ and then apply $(*)$ with $\mathcal A$ in place of $\mathcal S$ to find the conditional probability that all bins in $\mathcal A$ are full. Multiply this by the chance all items originally were in $\mathcal A.$

To illustrate, here is an R program to compute $(*).$ It uses sum, sum, prod, and rowSums to calculate the sum of sums of products of sums described by the foregoing formulas.

# A distribution on 1, 2, ..., n is represented by a vector of non-negative
# values summing to unity.  A set of distributions is a matrix whose rows
# represent distributions.
#
# Compute the chance all bins are full.
#
all.full <- function(p) {
N <- ncol(p)
sum(sapply(0:(N-1), function(n) {
(-1)^n * sum(apply(combn(1:N, N-n), 2, function(a) {
prod(rowSums(p[, a, drop=FALSE]))
}))
}))
}


Here is an example of its use.

p <- rbind(c(2,2,2), c(1,2,3), c(2,3,1), c(3,1,2))
p <- p / rowSums(p)
all.full(p)


The first two lines create an array whose $k=4$ rows represent discrete distributions on a set $\mathcal S$ of three elements. The output is $0.472222,$ a decimal representation of $17/36.$

As a check, here is a simulation of 50,000 iterations. The output gives the estimate and its standard error.

n <- ncol(p)
set.seed(17)
Sim <- (function(x) c(Mean=mean(x), SE=sd(x) / sqrt(length(x))))(replicate(5e4, {
length(setdiff(1:n, apply(p, 1, function(p) sample.int(n, 1, prob=p)))) == 0
}))


The result is

       Mean          SE
0.470020000 0.002232067


This is within one standard error of the result of all.full, suggesting the calculation and the simulation are in agreement.

• Thanks! How long do the 50k simulations run in R? – ondrejba Sep 11 '18 at 18:44
• A couple of seconds--if I post any code that is longer than that I warn my readers. It scales approximately linearly in $n$ and $k.$ – whuber Sep 11 '18 at 19:14