# Question

Consider a population of known size $N$, from which we sample $n$ individuals without replacement and measure their trait value $X_i$. All traits values are bounded between $[0,1]$. Let $\mu = \frac{\sum_i^N x_i}{N}$ be the average trait value in the population. If it makes things easier, I am happy to consider that $X$ is boolean ($0$ or $1$).

What is the unbiased estimator for $\mu - \mu^2$?

# Attempt

\begin{align} \mu - \mu^2 &= \operatorname{E}[\bar X] - \operatorname{E}[\bar X]^2\\ &= \operatorname{E}[\bar X] - \operatorname{E}[\bar X]^2 - \operatorname{var}[\bar X] + \operatorname{var}[\bar X] \\ &= \operatorname{E}[\bar X] - \operatorname{E}[\bar X]^2 - \operatorname{var}[\bar X] + \frac{1}{n} \sigma^2 \\ &= \operatorname{E}[\bar X] - \operatorname{E}[\bar X^2] + \frac{1}{n} \operatorname{E}[\hat \sigma^2], \end{align}

where $\hat \sigma^2$ is the estimator for the variance. Looking at this post,

$$\hat\sigma^2 = \frac{N-1}{N(n-1)} \sum_{i=1}^n (X_i-\bar X)^2$$

Hence,

$$\mu - \mu^2 = \operatorname{E}\left[ \bar X - \bar X^2 + \frac{1}{n} \frac{N-1}{N(n-1)} \sum_{i=1}^n (X_i-\bar X)^2\right]$$

# Test

Looked good to me so I tried it out (in R)

nbtrials = 5000      # number of independent samples
N = 20               # Number of individuals in the population
pop = rep(0:1, N/2)  # Create a boolean population with $\mu = 0.5$
n = 17               # Sample size

out = numeric(nbtrials)  # out will collect the estimates of $\mu - \mu^2$
for (trial in 1:nbtrials)
{
s = sample(pop,size=n, replace=FALSE)   # Take a sample without replacement
xbar = sum(s) / n                       # Compute the average frequency in the sample
out[trial] = xbar - xbar^2 + 1/n * (N - 1) / (N*(n-1)) * sum((s - mean(s))^2)   # Compute the estimator
}

# Now compare the population with the average of estimations taken from the samples
trueMean=sum(pop) / length(pop)
print(paste("True value  = ",trueMean *(1-trueMean)))
print(paste("Average estimated value = ",mean(out)))


which outputs something like

"True value  =  0.25"
"Average estimated value =  0.262343771626298"


showing a clear systematic overestimation. Where did I go wrong?

# Where does the mistake lie?

From @EinarRødland's comment, it seems likely that my mistake is in the equation $var[\bar X] = \frac{\sigma^2}{n}$, which is wrong due to the fact that sampling occurs without replacement from a finite population.

• sample(pop,size=n, replace=FALSE) I am suspicious of the foregoing command. 2000 times I iterated the command sum( sample(pop,size=n, replace=FALSE) ) and put all 2000 numbers in one long vector. The mean was $8.5225$; the standard deviation was $0.822$; the max and min were $10$ and $7$. What is this command doing? I don't see that it does anything like what you could have intended. – Michael Hardy Aug 12 '17 at 23:29
• It samples n individuals without replacement from the vector pop which is above build as containing only boolean values 0 and 1, each at frequency 0.5. The mean you got is around $n/2 = 8.5$ as expected. You could divide your sums by $n$ and the mean of the outputs should be close to mean(pop). I think it is doing what I think it is doing. Does it make sense to you? – Remi.b Aug 12 '17 at 23:55
• How does it do that? It gives a list of length $17$ in which each number is $0$ or $1$. How is that a choice of $17$ members of the list of $20\text{ ?}$ And why would the probability of being chosen by $0.5\text{?}$ Shouldn't you be trying to choose $17$ out of $20,$ rather than some random number out of $17\text{?}$ Moreover, instead of doing $5000$ random independent trials, couldn't you systematically list all $1140$ ways of choosing $17$ out of $20\text{?} \qquad$ – Michael Hardy Aug 13 '17 at 5:18
• I am not sure what is unclear for you about the function sample. In the code, it samples without replacement 17 elements from a vector of 20 elements. As, in the vector pop half of the elements are 1 (it is the arbitrary frequency of 1 that I chose in the population), the first sampled element has equal probability of being a zero or a one. You are right that it would have been smarter to consider the 1140 ways instead of making 5000 random trials. By taking 5000 random trials, it has the advantage to work even for when the number of ways is very large depending on the input parameters. – Remi.b Aug 13 '17 at 16:08
• This does NOT appear to be sampling without replacement 17 elements from a set of 20 elements. It gives me something that looks like this:  0 1 1 0 0 1 1 0 0 1 0 1 1 1 1 1 0. There are only 17 numbers, not 20. If it were sampling 17 out of 20, there should be 17 ones and 3 zeros, and the only randomness would be in choosing which three are zeros, and the probability $0.5$ would not be in any way involved. And your vector with 20 entries in which 0 and 1 simply alternate appears to have nothing to do with the problem. – Michael Hardy Aug 13 '17 at 18:29

$$\operatorname{var}[\bar X] = \frac{s^2}{n} \left(1-\frac {n-1} {N-1} \right) = \frac{\sum (X_i - \bar X)^2}{n(n-1)} \left(1-\frac {n-1} {N-1} \right)$$
$$\mu - \mu^2 = \operatorname E\left[\bar X - \bar X^2 + \frac{\sum (X_i - \bar X)^2}{n(n-1)} \left(1-\frac {n-1} {N-1} \right)\right]$$