# How to calculate $Var(X)$ for the union size?

Let's say there is a set of $$4$$ bags $$\{{a,b,c,d\}}$$ containing balls of colors $$\{{red,blue,green,orange,black\}}$$. Balls are assigned to bags in an arrangement that allows a variable number of colors in each bag, and overlaps.

$$a: \{{red, green\}}, b: \{{red, blue\}}, c: \{{orange\}}, d: \{{red, black\}}$$

When sampling $$k$$ bags, without replacement, one counts the total union of unique colors. Two examples:

• If $$k = 1$$, then one can sample either $$a$$, $$b$$, $$c$$, or $$d$$, obtaining $$\{{red, green\}}$$, $$\{{red, blue\}}$$, $$\{{orange\}}$$ or $$\{{red, black\}}$$ as colors, respectively. Then the number of unique colors is $$\{{red, blue\}}$$ = $$2$$ if selecting $$a$$, or $$\{{orange\}}$$ = $$1$$ if selecting $$b$$.
• If $$k = 2$$, one can draw $$a$$ and $$b$$, obtaining $$\{{red, blue\}}$$ and $$\{{red, blue\}}$$ as colors. Then the union size is the number of unique colors, $$\{{red, blue, orange\}}$$ = $$3$$.
1. What's $$Var(X)$$, when drawing $$k$$ bags from this set? Checking all options, one can calculate $$E(X) = 1.75$$ when $$k = 1$$, $$E(X) = 3$$ when $$k = 2$$ $$E(X) = 4$$ when $$k = 3$$, and $$E(X) = 5$$ when $$k = 4$$.
2. Is there a formulation for the underlying distribution that one can use, inferred from $$E(X)$$ or derived from all provided information? I think this could be a variant of a hypergeometric distribution, with variable weights per member, and duplicates.
• This would be easier with more concrete terminology. Eg: Suppose $n$ bags each have balls of one or more colors. Consider the total number of colors you get after randomly selecting $k$ bags. What is the variance of that total? Or: Suppose you have $n$ friends, each with some dietary restrictions. When you invite $k$ friends to an event with food, what is the variance of the total number of dietary restrictions you have to consider? Commented Jun 23, 2021 at 20:26
• Updated main body terminology (bags and colors). Commented Jun 23, 2021 at 20:43

The key to this solution is to represent the "bags" with a logical array (of zeros and ones) representing non inclusion of elements in the bags. Each bag gets its own row and each element of the universal set $$S = \{\text{black}, \ldots, \text{red}\}$$ gets its own column. This matrix will represent the "is not an element of" relation. That is, for a bag $$i$$ and element $$j$$, set $$a_{ij}=1$$ when $$j\notin i$$ and otherwise set $$a_{ij}=0$$ (when $$j\in i$$). Here is the matrix $$\mathbb{A}=(a_{ij})$$ for the example:

$$\mathbb{A} = \begin{array}{l|ccccc} \text{} & \text{black} & \text{blue} & \text{green} & \text{orange} & \text{red}\\ \hline \text{a} & 1 & 1 & 0 & 1 & 0 \\ \text{b} & 1 & 0 & 1 & 1 & 0 \\ \text{c} & 1 & 1 & 1 & 0 & 1 \\ \text{d} & 0 & 1 & 1 & 1 & 1 \end{array}$$

The sample space is the collection of all subsets of the bags $$\{a,b,c,d\}.$$ (It is okay for two or more bags to be identical as sets: everything still works out.)

We are going to define some random variables with values in $$\{0,1\}.$$ To do this, pick a column $$j$$ in the matrix and let the value of $$X_j$$ for any collection of $$\mathcal{B}=\{i_1,i_2,\ldots,i_k\}$$ be the product of all the entries for those bags:

$$X_j(\mathcal{B}) = \prod_{i\in\mathcal{B}} a_{ij}.$$

For instance, for the collection of bags $$\mathcal{B}=\{a,b,c\}$$ we will consider the corresponding rows of $$\mathcal{A}$$ (the first three). There are five corresponding random variables:

\begin{aligned} X_{\text{black}}(\{a,b,c\}) &= a_{a,\text{black}}\,a_{b,\text{black}}\,a_{c,\text{black}} = (1)(1)(1) = 1; \\ X_{\text{blue}}(\{a,b,c\}) &= a_{a,\text{blue}}\,a_{b,\text{blue}}\,a_{c,\text{blue}} = (1)(0)(1) = 0; \end{aligned}

and so on.

Clearly, the value of $$X_j$$ is $$0$$ if and only if $$j$$ is not an element of at least one bag in the collection. Consequently, letting $$n = |S|$$ be the number of columns of $$\mathbb{A},$$ the cardinality of the union of all bags in a collection $$\mathcal{B}$$ (its "union count") is given by

$$\text{Union count}(\mathcal{B}) = n - X(\mathcal{B})$$

where

$$X(\mathcal{B}) = \sum_{j\in S} X_j(\mathcal{B}).$$

I will analyze $$X,$$ because (this should be obvious) it will have the same variance as the union count.

On to the calculations. We need to find means, variances, and covariances of $$X$$ when the domain of $$X$$ is restricted to all subsets of size $$k.$$ (I will use a superscript $$(k)$$ to designate this restriction.) But this is now straightforward.

As a preliminary, for any column $$j$$ of $$\mathbb A$$ let $$m_j$$ be the sum of the column and for any pair of columns $$j, j^\prime$$ let $$m_{jj^\prime}$$ be the sum of the (componentwise) product of those columns:

$$m_j = \sum_{i} a_{ij};\quad m_{jj^\prime} = \sum_{i} a_{ij}\,a_{ij^\prime}.\tag{1}$$

(The latter are the entries of $$\mathbb{A}^\prime \mathbb{A}.$$)

Let $$m$$ be the number of rows of $$\mathbb A$$ (the number of bags).

• The mean of $$X_j^{(k)}$$ is the proportion of $$k$$-subsets of column $$j$$ of $$\mathbb A$$ where all the values are $$1.$$ Since there are $$m_j$$ ones in that column out of all $$m$$ entries, this proportion is the number of $$k$$-subsets of those ones, divided by the total number of possible subsets, $$E\left[X_j^{(k)}\right] = {\binom{m_j}{k}}/{\binom{m}{k}}.$$

• Because the values of every $$X_j^{(k)}$$ are in $$\{0,1\},$$ $$\left(X_j^{(k)}\right)^2 = X_j^{(k)}.$$ Thus, $$E\left[\left(X_j^{(k)}\right)^2\right] = {\binom{m_j}{k}}/{\binom{m}{k}}.$$

• Similarly, the expectation of $$X_j^{(k)}X_{j^\prime}^{(k)}$$ is found by counting all bags $$i$$ where both columns of $$\mathbb A$$ contain $$1$$'s. Thus $$E\left[X_j^{(k)}X_{j^\prime}^{(k)}\right] = \binom{m_{ij}}{k}/{\binom{m}{k}}.$$

We now have all the moments needed to compute the variance-covariance matrix

\begin{aligned} \mathbb{V^{(k)}}_{jj^\prime} &= \operatorname{Cov}\left(X_j^{(k)}, X_{j^\prime}^{(k)}\right) \\ &= \binom{m_{jj^\prime}}{k}/{\binom{m}{k}} - \left(\binom{m_{j}}{k}/{\binom{m}{k}}\right)\left(\binom{m_{j^\prime}}{k}/{\binom{m}{k}}\right)\\ &=\binom{m}{k}^{-2}\left(\binom{m}{k}\,\binom{m_{jj^\prime}}{k} - \binom{m_{j}}{k}\,\binom{m_{j^\prime}}{k}\right). \end{aligned}\tag{2}

Finally, because $$X^{(k)}$$ is the sum of the $$X_j,$$ its variance is

$$\operatorname{Var}\left(X^{(k)}\right) = (1,1,\ldots,1)\,\mathbb{V^{(k)}}\,(1,1,\ldots,1)^\prime = \sum_{j, j^\prime\in S}\mathbb{V^{(k)}}_{jj^\prime}.$$

This the sum of all the elements of this matrix, computed using the simple formulas $$(1)$$ and $$(2)$$ above.

I doubt there is much, if any, algebraic simplification available, because the variance must depend on the detailed relationship among all the bags: nothing short of making all bag-to-bag comparisons will do the trick.

The computational effort for this approach is proportional to the number of entries in the covariance matrix, on the order of $$n^2,$$ regardless of how many bags there might be. Furthermore, $$S$$ can be restricted to just the union of elements actually appearing in all the bags, potentially decreasing $$n.$$

The variances in the example are all zero except when $$k=1,$$ where the variance is $$3/16 = 0.1875.$$

The following R code implements this algorithm. After creating $$\mathbb A,$$ the calculations are noteworthy for their brevity. Functions mu1, mu2, and variance comprise just five lines of code. Using them, the calculation of the variances for all $$k$$ is simmply a matter of computing the variance matrix and summing its entries:

k <- 0:length(Bags)
V <- sapply(variance(k, A), sum)


Here's the full code. It includes a commented-out section to generate synthetic datasets and at the end it checks the calculation with a direct computation of the random variables $$X^{(k)}$$ and using the R function var to obtain the answers. I have run it with examples up to $$m=500$$ bags, where the calculation still takes less than a second. The brute-force method used to check the calculation is exponentially slow, however, and so is limited to small values of $$m.$$

#
# All first moments.
#
mu1 <- function(k, A) {
k <- pmax(0, pmin(k, nrow(A)))
exp(lchoose(colSums(A), k) - lchoose(nrow(A), k))
}
#
# All second moments.
#
mu2 <- function(k, A) {
lapply(k, function(k) matrix(exp(lchoose(crossprod(A), k) - lchoose(nrow(A), k)), ncol(A)))
}
#
# Variance matrix.
#
variance <- function(k, A) {
M <- lapply(k, function(k) (function(x) outer(x, x))(mu1(k, A)))
mapply(function(m2, m) m2 - m, mu2(k, A), M, SIMPLIFY = FALSE)
}
#------------------------------------------------------------------------------#
#
# EXAMPLES
#
# Represent the element relation for the example.
#
Bags <- list(a = c("red", "green"),
b = c("red", "blue"),
c = c("orange"),
d = c("red", "black"))
S <- sort(unique(unlist(Bags)))
# #
# # Alternatively, create random examples.
# # The algorithms resource requirements scale quadratically in the size of S,
# # which can be as large as n (but might be smaller).
# # It scales linearly with m, though!
# #
# n <- 5     # Number of elements in the universal set
# m <- 6    # Number of bags
# p <- .4     # Expected proportion in each bag
# S <- paste0("_", seq_len(n))
# # set.seed(17)
# Bags <- lapply(seq_len(m), function(i) S[runif(n) <= p])
# # Supply bag ids for up to 26^2+26=702 bags:
# s <- apply(expand.grid(c(letters), c("", letters)), 1, paste0, collapse="")
# names(Bags) <- s[seq_len(m)] #
#------------------------------------------------------------------------------#
#
# CREATE THE MATRIX A
#
S <- sort(unique(unlist(Bags)))
indices <- seq_along(S)
if (length(intersect(indices, S)) > 0) warning("Integer names of bags may confuse R.")
names(indices) <- S
one <- rep(1, length(S))
A.raw <- sapply(Bags, function(bag) {
a <- c(one, NA)        # Guarantees a non-empty result
a[indices[bag]] <- 0
a
})
A <- matrix(t(A.raw), ncol=length(S)+1)[,-(length(S)+1)] # Strips the NAs
colnames(A) <- S
rownames(A) <- names(Bags)
#------------------------------------------------------------------------------#
#
# The variances of the union counts.
#
k <- 0:length(Bags)
V <- sapply(variance(k, A), sum)
names(V) <- k
#------------------------------------------------------------------------------#
#
# CONFIRMATION
#
# Check with brute force.
#
if (length(Bags) <= 16) {
UnionCount <- function(a.list) length(unique(unlist(a.list)))
X <- sapply(k, function(k) {
i <- combn(rownames(A), k) # Columns are samples of bags w/o replacement
apply(i, 2, function(i) UnionCount(Bags[i]))
})
V.check <- sapply(X, function(x) {
n <- length(x)
ifelse(n <= 1, 0, var(x) * (n-1) / n)
})
names(V.check) <- k
#
# Compare.
#
print(results <- rbind(Formula = zapsmall(V), Brute force = V.check))
print(Bags)
print(A)
if(!isTRUE(all.equal.numeric(V, V.check))) stop("Discrepant results found.")
} else {
print(V)
print(A)
}

• This is great! I have tested it in R, and also a Python adaptation, and it runs very well versus permutations for $m < 500$. Tried approximating $nCk$ with Stirling's formula but still the main calculation bottleneck is $crossprod$. My use cases have $m$ values between $100$ and $1,000,000$, so it can be still used in a subset of cases. Thanks a lot! Commented Jun 25, 2021 at 16:58
• How many bags do you have? If it's small, the brute-force method can succeed. Otherwise you will need to resort to approximations.
– whuber
Commented Jun 25, 2021 at 17:02
• The number of bags varies between 100 and 100,000. The number of colors is also variable, and it can go between 10 and 1,000. Commented Jun 25, 2021 at 20:56
• Sampling the distribution for $m$ much greater than $1000$ or so might be the way to go.
– whuber
Commented Jun 25, 2021 at 21:30