# Estimating population size of a subgroup based on independent samples without replacement

Let a bag have 1000 balls, 100 red and 900 blue. Now, let us get ten independent samples (without replacement within each sample) of the original population (10% each). Proceed to count the number of red balls in each sample, i.e.:

(9, 10, 10, 11, 13, 8, 5, 15, 12, 9)


To estimate the number of red balls in the original population, one can calculate the average of the above list and divide by the sample. In this case:

(9 + 10 + 10 + 11 + 13 + 8 + 5 + 15 + 12 + 9) / 10 / 0.1 = 102


My questions are thus:

1. Instead of a number, how can I obtain the probability mass function of the estimated population? Is it sound to use the central limit theorem here, where the distribution $N = (\mu,\sigma^2/n)$ ? I have a problem with this approach since negative values don't make sense.
2. Now imagine I have an infinite numbers of samples (all w/ 10% of the population), or at least a very high (but known) number of samples, and keep getting the result for each of them, one by one. How do I decide when it is likely ok to stop (within a certain confidence)?

Some clarifications over the original problem:

1. The population size is known à priori.
2. But the number of colors is not known à priori.
3. Neither their ratio.
4. Each sample take balls without replacement.
5. Samples are independent (after a sample is taken, all balls are replaced).
• I don't have solid answers, but I can give pointers: for Question 1, you'd want to look into the binomial distribution (e.g., a rough 95% confidence interval in R: qbinom(p = c(0.025, 0.975), size = 1000, prob = mean(x)/100), where x is your sampled results). For 2), you're basically looking for a stopping rule; the clinical trials literature contains a lot of discussion on stopping rules. Commented Nov 21, 2013 at 20:04
• For Question 1, you should use a Bayesian approach with a beta prior for the proportion of red balls. For question 2, as you update your prior you will see the variance shrink as you continue to sample and you can stop based on whatever criteria you are comfortable with. Commented Nov 21, 2013 at 23:09
• Can you please clarify how you are performing the replacement? Do you mean, on each of 10 draws, that you are withdrawing 100 balls at a time (without individual replacement) and then replacing the entire 100 balls all together at the same time? Or are you drawing 1 ball individually and replacing it 100 times, and then looping over that process another 10 times? I suspect you mean the former, in which case your inner loop (over 100 balls) is technically taking place without replacement. If that's the case, the probability distribution you want is actually hypergeometric, not binomial. Commented Nov 23, 2013 at 7:10
• @stachyra: you are absolutely correct. Within each sample the balls are being taken without replacement, so in fact I might just be making 10 independent samples of the same population. Commented Nov 23, 2013 at 12:34
• @soakley: could you please clarify a little bit further? I don't know, nor believe in, any specific ratio of the balls à priori. Commented Nov 23, 2013 at 12:36

A Bayesian approach will work as follows:

1. Chose a prior distribution on the number $n$ of balls in the urn. If you don't have any prior idea, you can chose the uniform distribution : $\def\P{\mathbb P}$ $$\P_0(n = i) = {1\over 1001} \ \ (i = 0, \dots, 1000).$$

2. For a sample of 100 balls, denote $X$ the number of red balls. For each sample with $X=x$ update the distribution by using Bayes rule: \begin{align} \P_k(n = i) &= \P_{k-1}(n = i | X = x )\\ & = \P_{k-1} (n=i) \times { \P( X= x | n = i ) \over \sum_j \P( X=x | n = j ) \P(n=j) }. \end{align} Here $k$ is an index for the samples (starting from 1).

The probability $\P(X=x|n=i)$ is given by the hypergeometric distribution: $$\P(X=x|n=i) = { {i \choose x} {1000 -i \choose 100 - x } \over {1000\choose 100} }.$$

You then have to chose a stopping criteria. You can for example stop when there is a value which has a probability > 95% (this might take a while).

Here is some R code to perform this.

# flat prior
P0 <- rep(1/1001,1001)

# updating function
update <- function(P, x)
{
# P(X=x | n = i)
i <- 0:1000
Q <- dhyper(x, i, 1000-i, 100)
# P(X=x)
Px <- sum(Q*P)
return( Q*P/Px )
}

X <- c(9, 10, 10, 11, 13, 8, 5, 15, 12, 9)
P <- P0;
for(x in X) P <- update(P,x)

plot(P, type="h", xlim=c(50,150), col="red", lwd=2)