How to guess the size of a set?

Question

Assume we have a set of unique words and draw a number $n$ of them using simple-random-sampling without replacement independently in each round. We have several rounds and try to guess the set size after each round. The words are drawn from the same set in each round, but are drawn without replacement in each round (all words are distinct in each round). How could this be achieved and how big would the error be?

For example if the set size is smaller than $n$, then we would only need one round. If the set size is close to $n$, in the second round the number of words that were drawn in the first round would be close to $n$. So I know that we get some information, but I don't understand how we could efficiently guess the set size.

Answer 1

First some notation: Assume $N$ unique words, sampled with the same probability. We are doing $R$ (independent) rounds of simple random sampling without replacement, each round the sample size is $n \le N$ (generalization to unequal sample sizes should be straightforward).

Define indicator random variables $$ X_{ij} =\begin{cases} 1 & \text{word $j$ sampled in round $i$} \\ 0 & \text{otherwise} \end{cases} $$ Note that we have \begin{align} \sum_j X_{ij} &= n ~(\text{for all $i$}) \\ \sum_i X_{ij} &\sim \mathcal{Binom}(R, n/N) \end{align} but all these binomial random variables are not independent, since they have a fixed, constant sum $nR$. But if $N$ is large and $n/N$ small the dependence would be slight.

Define $Y_j = \sum_i X_{ij}$ each having the binomial distribution defined above. The exact likelihood function for this problem will be intractable, so I will not write it out. But based on this binomials $Y_j$ we can construct a composite likelihood function, see the references at Parameter Estimation for intractable Likelihoods / Alternatives to approximate Bayesian computation.

But before going into the details, we can also define pair statistics $$ Y_{jl} = \sum_{i=1}^R X_{ij} X_{il} \quad \text{for $j\not=l$} $$ which will also have binomial distributions: $\mathcal{Binom}(R,\frac{n}{N}\cdot\frac{n-1}{N-1})$.

The idea with composite likelihood (a specific form of pseudo-likelihood) is to construct individual likelihood functions from parts of the data, and then just multiply them together, as would be correct if they where independent ... even if they are not independent. See the linked page above for details.

First, we find the composite likelihood only based on the $Y_j$. First, suppose the total number of unique words sampled in the $R$ rounds are $M \le N$. Then the composite likelihood for the one unknown parameter $N$ becomes $$ \prod_{j=1}^M \binom{R}{y_j}(n/N)^{y_j} (1-n/N)^{R-y_j} \cdot \left[ (1-n/N)^R \right]^{(N-M)} $$ the last factor coming from the $N-M$ unobserved words.

(sorry, now it is late night here, so I will continue this answer tomorrow)

Answer 2

Using the following notation:

$s_i$: the number of samples for the $i^{th}$ round

$k_i$: the number of words sampled in the $i^{th}$ round that had not been previously sampled

$m_i=\sum_{j=1}^ik_j$

$\pi_1(n)$: the prior distribution of $n$

We'll also assume $n\geq{s}$; otherwise we can determine $n$ after the first round of sampling. We can update $\pi(n)$ beginning with the second round of sampling:

$$\pi_i(n)\propto\frac{\binom{m_{i-1}}{s_i-k_i}\binom{n-m_{i-1}}{k_i}}{\binom{n}{s_i}}\pi_{i-1}(n)$$ $$\propto\pi_1(n)\prod_{j=2}^i\frac{(n-m_{j-1})!(n-s_j)!}{n!(n-m_{j-1}-k_j)!}$$

Implemented as an R function:

pn <- function(s, k, prior) {
  l <- length(s)
  m <- cumsum(k[-l])
  s <- s[-1]
  k <- k[-1]
  function(n) pmax(prior(n)*exp(colSums(lgamma(outer(1 - m, n, "+")) + lgamma(outer(1 - s, n, "+")) - lgamma(outer(1 - m - k, n, "+"))) - (l - 1)*lgamma(n + 1)), 0, na.rm = TRUE)
}

For example, say $\pi_1(n)\sim{U(8,30)}$, $k=5,3,1$, and $s=5,5,5$.

k <- c(5, 3, 1)
s <- rep(5, 3)

post <- pn(s, k, function(n) 1)
like <- post(8:30)
plot(8:30, like/sum(like), xlab = "n", ylab = "p(n)", col = "blue", pch = 3)

We can verify the results with simulation:

library(parallel)
set.seed(724526144)
nreps <- tabulate(sample(23, 23e6, TRUE), 23)
clust <- makeCluster(detectCores() - 1)
clusterExport(clust, c("nreps", "s", "k"))
sim <- unlist(parLapply(clust, 1:23, function(i) sum(replicate(nreps[i], all(cumsum(!duplicated(c(replicate(3, sample(i + 7, 5)))))[cumsum(s)] == cumsum(k))))))
stopCluster(clust)
points(8:30, sim/sum(sim), col = "orange")
legend("topright", legend = c("Posterior likelihood", "Simulation"), col = c("blue", "orange"), pch = c(3, 1))

If we had observed $k=5,4,2$ instead of $k=5,3,1$, we can see the posterior distribution shift toward larger values of $n$.

k <- c(5, 4, 2)

post <- pn(s, k, function(n) 1)
like <- post(8:30)
plot(8:30, like/sum(like), xlab = "n", ylab = "p(n)", col = "blue", pch = 3)

Again verifying the results with simulation:

set.seed(353678169)
nreps <- tabulate(sample(23, 23e6, TRUE), 23)
clust <- makeCluster(detectCores() - 1)
clusterExport(clust, c("nreps", "s", "k"))
sim <- unlist(parLapply(clust, 1:23, function(i) sum(replicate(nreps[i], all(cumsum(!duplicated(c(replicate(3, sample(i + 7, 5)))))[cumsum(s)] == cumsum(k))))))
stopCluster(clust)
points(8:30, sim/sum(sim), col = "orange")
legend("topright", legend = c("Posterior likelihood", "Simulation"), col = c("blue", "orange"), pch = c(3, 1))

Answer 3

This is known as the Mark and Recapture problem, because the people who most often encounter it are ecologists estimating population size. There are many methods available, which you can find by searching for "Mark and Recapture".