2
$\begingroup$

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.
$\endgroup$
2
  • $\begingroup$ 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? $\endgroup$
    – Matt F.
    Commented Jun 23, 2021 at 20:26
  • $\begingroup$ Updated main body terminology (bags and colors). $\endgroup$
    – ilibarra
    Commented Jun 23, 2021 at 20:43

1 Answer 1

1
$\begingroup$

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)
}
$\endgroup$
4
  • 1
    $\begingroup$ 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! $\endgroup$
    – ilibarra
    Commented Jun 25, 2021 at 16:58
  • $\begingroup$ How many bags do you have? If it's small, the brute-force method can succeed. Otherwise you will need to resort to approximations. $\endgroup$
    – whuber
    Commented Jun 25, 2021 at 17:02
  • $\begingroup$ 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. $\endgroup$
    – ilibarra
    Commented Jun 25, 2021 at 20:56
  • 1
    $\begingroup$ Sampling the distribution for $m$ much greater than $1000$ or so might be the way to go. $\endgroup$
    – whuber
    Commented Jun 25, 2021 at 21:30

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.