# To find variance and covariance for a double sampling problem

A simple random sample of size $n=n_1 + n_2$ is drawn without replacement from a finite population of size $N$. Further a simple random sample of size $n_1$ is drawn without replacement from the first sample. Let $\bar y$ and $\bar {y_1}$ be the respective sample means.

Find V($\bar {y_1}$) and V ($\bar {y_2}$), where $\bar {y_2}$ is the mean of the remaining $n_2$ units in the first sample. Also find Cov($\bar {y_1}$,$\bar {y_2}$). Assume that $S^2$ is the population variance.

This is what I've been able to do thus far:

Note that $\bar {y_2}$ = $\frac{(n\bar y) - n_1\bar {y_1}}{n_2}$.

Also, V($\bar {y_1}$) = $E_1[V_2(\bar {y_1})] + V_1[E_2(\bar {y_1})]$. After simplification, this has turned out to be $(\frac1{n_1} - \frac1N)S^2$. Now am I on the right track? Also, some help in finding $V(\bar y_2)$ and the covariance would be appreciated.

• Hint: you already know $\text{V}(\bar{y}_2)$: the formula is in your question.
– whuber
May 7, 2014 at 15:14
• @whuber, kindly point it out meticulously. It is turning out to be quite complicated. May 7, 2014 at 17:40
• Interchanging $1$ and $2$ in your question produces the answer for $\text{V}(\bar{y}_2)$, because the complement of a random sample without replacement is itself a random sample without replacement. This leaves the hard part, which is finding the covariance.
– whuber
May 7, 2014 at 17:54

Following the analysis of a single sample without replacement at https://stats.stackexchange.com/a/622287/919, the variance of the mean of any sample of size $$K$$ from the population $$(x_1, x_2, \ldots, x_N)$$ with variance $$\sigma^2$$ (called $$S^2$$ in the present question) is

$$\operatorname{Var}(\bar X) = \frac{N-K}{K(N-1)}\,\sigma^2.$$

In the present double sampling scenario, for each subject index $$i$$ let $$J_i$$ indicate whether $$x_i$$ is in the second sample of size $$n_1$$ and $$I_i$$ indicate whether $$x_i$$ is in the remainder of the sample. Translating the description in the question into properties of these random variables gives

$$n_1 = \sum_{i=1}^N J_i,\tag{1a}$$ $$n_2 = \sum_{i=1}^N I_i,\tag{1b}$$ and $$I_iJ_i = 0,\ i = 1, 2, \ldots, N.\tag{1c}$$

We need to compute the second moments (variances and covariances) of these $$2N$$ variables. As noted in the previous thread, the $$(I_i,J_i)$$ are exchangeable--all pairs of nested subsets of the population are equally likely to result from this double sample--and we can exploit this in the calculation.

We already know (as shown in the initially referenced post) that $$(1a)$$ and $$(1b)$$ imply that for $$i\ne j,$$

$$\operatorname{Cov}(J_i,J_j) = -\frac{n_1(N-n_1)}{N^2(N-1)}$$

with a similar expression for $$\operatorname{Cov}(I_i,I_j).$$ The only new calculations needed are for the cross-covariances of $$I_i$$ and $$J_j.$$ Consider the expectation of the product of $$(1a)$$ and $$(1b)$$ and use $$(1c)$$ to simplify the diagonal terms:

\begin{aligned} n_1n_2 &= E\left[\sum_{j=1}^{N} J_j \sum_{i=1}^N I_i\right]\\ &= E\left[\sum_{i=1}^{N} I_iJ_i\right] + E\left[\sum_{i\ne j}^N I_iJ_j\right]\\ &= \sum_{i=1}^N E[0] + \sum_{i\ne j}^N E[I_iJ_j]\\ &= N(N-1)E[I_1J_2] \end{aligned}

The last equality is consequence of the exchangeability.

Solving, we find that for all $$i\ne j,$$

$$E[I_iJ_j] = \frac{n_1n_2}{N(N-1)}.$$

A standard formula for covariance shows (still for $$i\ne j$$)

$$\operatorname{Cov}(I_i, J_j) = E[I_iJ_j] - E[I_i]E[J_j] = \frac{n_1n_2}{N(N-1)} - \frac{n_1}{N}\frac{n_2}{N} = \frac{n_1n_2}{N^2(N-1)}.$$

The remaining calculation that is needed is simple:

$$\operatorname{Cov}(I_i,J_i) = E[I_iJ_i] - E[I_i]E[J_i] = -\frac{n_1n_2}{N^2}.$$

Now we can proceed without further thought--entirely mechanically--to find the covariance of the subsample means. Let $$\bar X = \frac{1}{n_2}\sum_{i=1}^N I_i x_i,\quad \bar Y = \frac{1}{n_1}\sum_{j=1}^N J_j x_j$$

be the means of the remainder and of the subsample. Let's clear the fractions and work directly with the sums, writing $$\mu$$ for the population mean (which, as a check of the preceding work, had better not enter into the result!):

\begin{aligned} n_1n_2\operatorname{Cov}(\bar X, \bar Y) &= \operatorname{Cov}(n_2\bar X, n_1\bar Y) = \operatorname{Cov}\left(\sum_{i=1}^N I_i x_i \sum_{j=1}^N J_j x_j\right)\\ &= \sum_{i=1}^N x_i^2 \operatorname{Cov}(I_i,J_i) + \sum_{i\ne j} x_ix_j \operatorname{Cov}(I_i,J_j)\\ &= -\frac{n_1n_2}{N^2}\sum_{i=1}^N x_i^2 + \frac{n_1n_2}{N^2(N-1)} \sum_{i\ne j}x_ix_j\\ &= \frac{n_1n_2}{N}\left[-\frac{1}{N}\sum_{i=1}^N x_i^2 + \frac{1}{N(N-1)}\left(\sum_{i,j=1}^N x_ix_j - \sum_{i=1}^N x_i^2\right)\right]\\ &= \frac{n_1n_2}{N}\left[-(\sigma^2 + \mu^2) + \frac{1}{N(N-1)}\left((N\mu)^2 - N(\sigma^2+\mu^2)\right)\right]\\ &= \frac{n_1n_2}{N}\left[\frac{-N}{N-1}\,\sigma^2\right]\\ &= -\frac{n_1n_2}{N-1}\,\sigma^2. \end{aligned}

Thus

$$\operatorname{Cov}(\bar X, \bar Y) = -\frac{\sigma^2}{N-1}.$$

Incidentally, we didn't need much of the preliminary algebra. We could have solved for $$\operatorname{Cov}(I_i,J_j),$$ where $$i\ne j,$$ as the unique value that results in the coefficient of $$\mu^2$$ equalling zero.

I was surprised that this result does not depend on the subsample sizes $$n_1$$ and $$n_2$$ or even on their sum, so I tested it with large amounts of simulated data. It looks correct: in the following R code, the covariance matrix estimated from 10,000 independent iterations of this experiment is very close to the variances found here. It outputs the term-by-term ratio of those two matrices, so at a glace (by comparing it to the matrix of $$1$$s) we can see how accurate the formulas are. Here's the output with the seed set to 17:

     [,1] [,2]
[1,] 1.02 1.02
[2,] 1.02 1.00


The chi-squared distribution tells us these ratios are likely to lie between 0.97 and 1.03, and so they are. Similar results ensue when the $$n_i$$ are varied (in the vector n), demonstrating the lack of dependence on the subsample sizes.

#
# Create a population.
#
N <- 20      # Population size
n <- c(5, 8) # Subsample sizes: sum cannot exceed N
# set.seed(17)

P <- sample.int(2*N, N, replace = TRUE)
sigma2 <- mean((P -  mean(P))^2)         # The population variance
#
# Sample and subsample.
#
n.iter <- 1e4
xy.bar <- replicate(n.iter, {
i <- sample.int(N, sum(n))       # Sample n[1] + n[2] indexes without replacement
j <- sample.int(length(i), n[1]) # Subsample n[1] indexes from the sample
y <- P[i[j]]                     # The subsample
x <- P[i[-j]]                    # The remainder
c(mean(y), mean(x))              # Their means
})
#
# Estimate the variance-covariance matrix relative to the population variance.
#
S <- var(t(xy.bar)) / sigma2
#
# Compare to the theoretical formula.
#
Sigma <- matrix(c((N - n[1])/(n[1] * (N-1)),
rep(-1 / (N-1), 2),                 # Doesn't depend on n!
(N - n[2])/(n[2] * (N-1))), 2, 2)

ndigits <- -log10((qchisq(1 - 0.05/2, n.iter - 1) / n.iter) - 1) # 95% confidence
(round(S / Sigma, ndigits))