# Estimating Size of a Set based on two Overlapping Subsets

I've searched everywhere for a similar question and many things come close but are not the same.

I'm looking for a way to estimate the size of a set if two partially overlapping subsets are known (assuming both subsets were selected at random).

Here's the simplified scenario:
One person randomly selects 100 numbers from the phone book and writes them down. A second person writes down 150 phone numbers. You obtain both lists, and find that 75 of the phone numbers are the same. Is there a way to estimate the total number of numbers in the phone book? Obviously there must be AT LEAST 175, but there are probably more.

I suppose you could imagine the same phrased as a lottery with non-sequentially numbered balls. After two drawings, you know a certain number of balls came up in both drawings, so can you estimate the total number of balls being drawn from?

I imagine there's some Bayesian solution, but I'm not sure.

• I took the liberty of adding a couple of relevant tags to help others find this question. Feb 12, 2015 at 1:05

This sounds like the basic "capture-recapture" problem, sometimes called "mark and recapture".

You have a population of unknown size $N$; imagine them to be indistinguishable balls in an urn (all white, say). You take a sample of size $n$ randomly from the population, and mark them (paint them black say), return them to the population, and mix.

You then draw a new sample, of size $m$, of which $k$ are marked.

This is of course a hypergeometric model (i.e. $k$ is hypergeometric).

The aim here is to estimate $N$.

(In your example, the set of numbers selected by the first person are the "marked" ones.)

You can use a variety of methods to estimate $N$.

Note that the mean of the hypergeometric is $mn/N$, so a naive method of moments estimate is $\hat{N}=mn/k$. In your example, you'd guess $N=200$. In the capture-recapture literature this is called the Lincoln–Petersen estimator. It's intuitively appealing because it equates sample proportion and population proportion; asymptotically, the sample proportion will converge to the population proportion.

Obviously, since in some cases $k$ can be 0, the estimator can (with non-zero probability) be non-finite, which is somewhat of a bias problem (indeed, $E(\frac{mn}{k})>N$ even if $k$ can't be zero); if you modify your estimator when $k$ is quite small, it can nevertheless perform fairly well).

An estimator that is notionally similar is the Chapman estimator $\hat{N} = \frac{(K+1)(n+1)}{k+1} - 1$. It performs substantially better in small samples.

Note that both of these estimators can yield noninteger estimates.

The maximum likelihood estimator: the likelihood is increasing* in $N$ for integers below $\frac{mn}{k}$ and increasing for integers above it; the integers $\lfloor \frac{mn}{k} \rfloor$ and $\lfloor \frac{mn}{k} \rfloor+1$ would seem to be the two possible candidates for maximizing the likelihood. It's a simple matter to directly compute the likelihood for both.

* in some circumstances it's actually nondecreasing between the last two points

Here's the likelihood function for a small interval around the method of moments estimator: It turns out the likelihood is equally high for $\hat{N}=199$ and $\hat{N}=200$. (Indeed, since the $mn/k$ is biased upward, it might make some sense to choose the lower of the two in this instance, $\hat{N}=199$.)

According to Zhang (2009), $\lfloor mn/k\rfloor$, the integer part of the method of moments estimator maximizes the likelihood (i.e. always round down). (I haven't checked this, but the argument looks sound. It doesn't hurt to directly compare a couple of values around $mn/k$ in any case -- in our example we discovered that the next lower value also maximizes the likelihood, though I think that can only happen when $mn/k$ is integer.)

Bayesian estimation is quite useful in this problem (when $k$ is small there's not a great deal of information in the sample to hold the "tail" down, so prior information about population size can be very useful), but of course the posterior distribution depends on the particular prior one chooses, and the estimator itself depends on the loss function selected.

Confidence intervals

Let's say we want a confidence interval for $N$, along with our estimate $\hat{N}$.

We can form a large-sample interval easily enough.

We know that $k$ is hypergeometric. From the normal approximation to the hypergeometric, $k$ is approximately $\sim N\left(\frac{mn}{N},\frac{mn(N-m)(N-n)}{N^2(N-1)}\right)\,$.

So $(k-\frac{mn}{N})/\sqrt{\frac{mn(N-m)(N-n)}{N^2(N-1)}}) = (kN-mn)/\sqrt{\frac{mn(N-m)(N-n)}{(N-1)}}$ is approximately standard normal.

By squaring and comparing with the upper 95% point of a $\chi^2_1$, we get an inequality in $N$ that can be rearranged into a cubic inequality: given $m,n$ and $k$, we find those values of $N$ such that $(kN-mn)^2(N-1)-3.84(mn(N-m)(N-n)) < 0$.

For the example problem that yields the following cubic: (Being a cubic, there's another interval where the function is negative, but that's in an impossible region for $N$ in this instance)

The cubic curve crosses the horizontal axis at about N=186.75 and N=219.78

This suggests that an approximate 95% interval for $N$ should be something like $(187,220)$, though you might "round outward" to be safer.

We could possibly derive an explicit solution for the approximate bounds, but unless you're doing many hundreds of these it's probably not worth the effort to do more than find the zeros fairly automatically (polynomial or even general root-finding functions are widely available).

In practice one should also investigate the actual coverage of intervals generated in this fashion for values of $N$ in the region of the estimate. That is, given $m$ and $n$, choose an $N$, simulate many $k$ (say 1000 or 10000), and find the proportion of intervals that include that $N$, repeating the exercise at several plausible values of $N$.

In the case of the Chapman estimator, we can obtain an approximate variance estimate for $\hat{N}$:

$\operatorname{var}(\hat{N}^{_\text{C}}) = \frac{(m+1)(n+1)(m-k)(n-k)}{(k+1)(k+1)(k+2)}$

(The $C$ superscript is to indicate the Chapman estimator)

In this case, $\hat{N}^{_\text{C}}$ should be asymptotically normal, so an asymptotic 95% interval would be

$\hat{N}^{_\text{C}}\pm 1.96 \sqrt{\frac{(m+1)(n+1)(m-k)(n-k)}{(k+1)(k+1)(k+2)}}$

which on your example gives $(184.95,216.4)$.

For the MLE, you'd get an asymptotic interval by approximating the log-likelihood at the peak by a quadratic (in effect, an asymptotic normal approximation for $\hat N$) to estimate $\hat{\sigma_{\hat N}}$ (which can be obtained from the approximating quadratic's second derivative) and from that normal approximation derive confidence limits. I won't labor the point, as details of such calculations can be found for a variety of MLEs, and are similar here.

: Zhang, H. (2009),
"A Note About Maximum Likelihood Estimator in Hypergeometric Distribution,"
Comunicaciones en Estadística, June, Vol. 2, No. 1

• Just to complicate things, what if the second person cheated a little... they actually wrote down 125, then copied down 25 numbers from the first person. I'm guessing this would be solved by only solving the same, but reducing both the overlap and person 2's count by 25... i.e. n=100, m=125, k=50; so thus N=250 by the simple hypergeometric calculation. Am I doing this properly? Feb 12, 2015 at 2:08
• I don't think that's the appropriate solution. You might want to ask that as a new question, though, since it moves away from the 'classic' problem addressed here. Feb 12, 2015 at 2:12
• Sure, I can ask that as a separate question. What's the correct way to state the uncertainty in this estimation? Feb 12, 2015 at 17:47
• I'm not sure what you're asking for there. Are you asking for some sort of uncertainty in the estimates of $N$ in the answer above? Feb 12, 2015 at 21:08
• What if there are 3 or more subsets, not just two? Jul 30, 2018 at 1:33