Capture-recapture method to estimate population size

caprecap.py

# Estimate the population size using capture-recapture method
# Assumes: experiment twice, closed population

import numpy as np
import matplotlib.pyplot as plt

T = 1000  # true pop size
N1 = 50  # size of sample 1
N2 = 30  # size of sample 2

population = np.arange(T)
ests = []
hit = 0

for i in range(100):

    # First sample + tag (the number it self) + mix back
    s1 = np.random.choice(population, size=N1, replace=False)
    # Second sample
    s2 = np.random.choice(population, size=N2, replace=False)

    # Recaptured
    R = len(set(s1).intersection(s2))
    # Estimate for pop size and variance of estimator
    if R > 0:
        mu = N1 * N2 / R
        var = (N1 * N2 * (N1 - R) * (N2 - R)) / R**3
    else:
        n1 = N1 + 1
        n2 = N2 + 1
        r = R + 1
        mu = n1 * n2 / r - 1
        var = (n1 * n2 * (n1 - r) * (n2 - r)) / r**3

    # 95% conf interval
    conf = 1.96 * np.sqrt(var)
    print(f"{mu:.2f} +/- {conf:.2f}")
    ests.append(mu)
    hit += abs(T - mu) <= conf

print(hit / 100)

plt.hist(ests)
plt.show()