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.
 A: 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.
