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)
}
