Finding the covariance of two random sums

Question

I am trying to derive the covariance of two sample sums.

Some notation and details:

$x_i$ and $y_i$ are numeric values of two characteristics corresponding to member i of a finite population of N members.

$w_i$ is a random indicator variable taking a value of 1 with a probability of n/N, and is otherwise 0. All of the N $w_i$ values sum to n.

So essentially, I’m taking a sample of size n from a population of size N, and I’m separately summing up the x and the y values for that sample. I want to find the covariance of these two sums, and I think I’m close but could use some help.

$$ Cov\left(\sum_{i=1}^{N}{w_ix_i},\sum_{i=1}^{N}{w_iy_i}\right) $$

---

*** Edit Following is my original post (with incorrect derivation) but scroll down for the correct answer including some R code confirming the answer. For a derivation of the answer, see the extremely helpful post by whuber below.

---

$$ =\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-E\left(\sum_{i=1}^{N}{w_ix_i}\right)\right)\left(\sum_{i=1}^{N}{w_iy_i}-E\left(\sum_{i=1}^{N}{w_iy_i}\right)\right)\right] $$ $$ =\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-\sum_{i=1}^{N}{x_iE\left(w_i\right)}\right)\left(\sum_{i=1}^{N}{w_iy_i}-\sum_{i=1}^{N}{y_iE\left(w_i\right)}\right)\right] $$ $$ =\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-\frac{n}{N}\sum_{i=1}^{N}x_i\right)\left(\sum_{i=1}^{N}{w_iy_i}-\frac{n}{N}\sum_{i=1}^{N}y_i\right)\right] $$$$ =\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-n\mu_x\right)\left(\sum_{i=1}^{N}{w_iy_i}-n\mu_y\right)\right]$$$$ =\ E\left[\sum_{i=1}^{N}{w_ix_i}\sum_{i=1}^{N}{w_iy_i}-\ n\mu_y\sum_{i=1}^{N}{w_ix_i}-\ n\mu_x\sum_{i=1}^{N}{w_iy_i}+n^2\mu_x\mu_y\right]$$$$ =\ E\left[\sum_{i=1}^{N}{w_ix_iw_iy_i}\right]+E\left[\sum_{i=1}^{N}\sum_{j\neq i}^{N}{w_ix_iw_jy_j}\right]-\ n\mu_y\frac{n}{N}\sum_{i=1}^{N}x_i-\ n\mu_x\frac{n}{N}\sum_{i=1}^{N}y_i+n^2\mu_x\mu_y$$$$ =\ \left(\sum_{i=1}^{N}{x_iy_i}\right)\left(\frac{n}{N}\right)+\left(\sum_{i=1}^{N}\sum_{j\neq i}^{N}{x_iy_j}\right)\left(\frac{n}{N}\right)\left(\frac{n-1}{N-1}\right)-\ n^2\mu_y\mu_x-\ n^2\mu_x\mu_y+n^2\mu_x\mu_y$$$$ =\ n\mu_{xy\left(i=j\right)}+n\left(n-1\right)\mu_{xy\left(i\neq j\right)}-\ n^2\mu_y\mu_x$$$$ =\ n^2\left(\frac{\mu_{xy\left(i=j\right)}+\left(n-1\right)\mu_{xy\left(i\neq j\right)}}{n}-\mu_x\mu_y\right)$$$$ =\ n^2(\mu_{xy}-\mu_x\mu_y)$$ $$=\ n^2Cov(x_i,y_i)$$

The above looks right to me, but I must be doing something wrong. When I simulate the covariance in R, I find that it is not the equation above. However, it seems (consistent with simulation) I can obtain the correct covariance by applying a finite population correction and by adding an n to the denominator of my result, as shown below.

$$\frac{N-n}{N}nCov(x_i,y_i)$$

---

*** End of original post: Following is the answer and some R code that "confirms" it.

---

$$\frac{n(N-n)}{N-1}Cov(x_i,y_i)$$ (where Cov() represents the population covariance)

R Code with results that are consistent with this answer:

library(mvtnorm)
library(tidyverse)

N <- 5
n <- 3

val1 <- rmvnorm(n = N, mean = c(50, 100), sigma = matrix(c(15^2, 12^2, 12^2, 15^2), nrow = 2))
colnames(val1) <- c("var1","var2")

set.seed(83442)

numSims <- 100000
covList1 <- as.list(1:numSims)
for(i in 1L:numSims) {
  covList1[[i]] <- val1[sample.int(N, n),] %>%
    as_tibble() %>%
    summarize(var1 = sum(var1), var2 = sum(var2))
}

covDF1 <- do.call("bind_rows", args = covList1)

(cov1 <- cov(covDF1$var1, covDF1$var2)) # Covariance from simulations: 195.7401
((N-n)/(N-1))*n*cov(val1[, 1], val1[, 2])*(N-1)/N # Covariance with N-1 denominator: 196.2536
# Note that the (N-1)/N at the end adjusts for the fact that we need a population (not sample) covariance. 

The results of this comparison improve with more simulations.

Answer 1

To help make the ideas clear, I will use capital letters for random variables.

Everything follows from the restriction $\sum W_i=n,$ because that implies this sum has zero variance. Since each $W_i$ is a Bernoulli variable,

$$\operatorname{Var}(W_i) = \frac{n(N-n)}{N^2}.$$

Computing the variance of the sum and assuming, as is the case with simple random sampling, that for $i\ne j$ $\operatorname{Cov}(W_i,W_j)$ does not depend on $i$ or $j,$ we find

$$\begin{aligned} 0 &= \operatorname{Var}\left(\sum_{i=1}^N W_i\right) \\ &=\sum_{i=1}^N \operatorname{Var}\left(W_i\right) + \sum_{i\ne j}^N \operatorname{Cov}(W_i,W_j) \\ &= N\frac{n(N-n)}{N^2} + N(N-1)\operatorname{Cov}\left(W_1,W_2\right), \end{aligned}$$

enabling us to solve for the covariance as

$$\operatorname{Cov}\left(W_i,W_j\right) = \operatorname{Cov}\left(W_1,W_2\right) = -\frac{n(N-n)}{N^2(N-1)}.$$

Consequently, assuming $N\gt 1,$ for fixed coefficients $(x_i)$ and $(y_i)$ and writing $\bar x = \sum x_i/N,$ $\bar y = \sum y_i/N,$ and $\overline{xy}=\sum_{i}x_iy_i/N,$ we find

$$\begin{aligned} \operatorname{Cov}\left(\sum_{i=1}^N x_iW_i, \sum_{j=1}^N y_jW_j\right) &= \sum_{i=1}^N x_iy_i \operatorname{Var}\left(W_i\right) + \sum_{i\ne j}^N x_iy_j\operatorname{Cov}\left(W_i,W_j\right) \\ &= \frac{n(N-n)}{N^2}\sum_{i=1}^Nx_iy_i - \frac{n(N-n)}{N^2(N-1)}\sum_{i\ne j}^N x_iy_j \\ &= \frac{n(N-n)}{N}\overline{xy} - \frac{n(N-n)}{N-1} \bar{x}\bar{y} + \frac{n(N-n)}{N(N-1)}\overline{xy}\\ &= \frac{n(N-n)}{N-1}\left(\overline{xy} - \bar{x}\bar{y} \right). \end{aligned}$$

(When $N=1$ the double sum doesn't appear and the result easily reduces to $0.$)

If we draw one of the $(x_i,y_i)$ randomly and equiprobably from all $N$ of these paired values, we have a bivariate random variable $(X,Y),$ enabling the result to be written

$$\operatorname{Cov}\left(\sum_{i=1}^N x_iW_i, \sum_{j=1}^N y_jW_j\right) = \frac{n(N-n)}{N-1} \operatorname{Cov}(X,Y).$$

---

I was tempted to check this result with simulation, but elected to use an exhaustive enumeration of all the possible samples instead, of which there are $\binom{N}{n}.$ For small $N$ this is feasible and gives precise results. The output computes the covariance of the weighted sums in three ways: using the formula in terms of $\overline{xy}-\bar{x}\bar{y},$ the formula in terms of $\operatorname{Cov}(X,Y),$ and--this is the verification--the population covariance of all possible sample sums.

An example of its output for $N=20,$ $n=15$ is

Direct Formula Covariance formula         Exhaustive 
     -3.035239          -3.035239          -3.035239 

demonstrating agreement in this case.

Here's the R code.

#
# Create *any* bivariate population you like.
#
N <- 20
# set.seed(17)
x <- rnorm(N)
y <- rnorm(N) - x
#
# Specify the sample size.
#
n <- 15
if(choose(N, n) > 1e6) stop("Are you sure you want to do this?", call.=FALSE)
# 
# Compute the distribution of the sample sum.
#
W <- combn(1:N, n)
wx <- apply(W, 2, function(w) sum(x[w]))
wy <- apply(W, 2, function(w) sum(y[w]))
#
# Compare various formulae.
#
c(`Direct Formula` = n * (N-n) / (N-1) * (mean(x*y) - mean(x)*mean(y)),
`Covariance formula`=n * (N-n) / N * cov(x, y),
Exhaustive = cov(wx, wy)*(1 - 1/length(wx)))

# plot(wx, wy) # Can be interesting...
```

Answer 2

Based on some of the methods whuber used in his answer, I decided to derive this covariance again, but this time in the way that I originally attempted the derivation (starting by showing the covariance as an expectation and going from there). I get the right answer now:

$$Cov\left(\sum_{i=1}^{N}{w_ix_i},\sum_{i=1}^{N}{w_iy_i}\right)$$ $$=\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-E\left(\sum_{i=1}^{N}{w_ix_i}\right)\right)\left(\sum_{i=1}^{N}{w_iy_i}-E\left(\sum_{i=1}^{N}{w_iy_i}\right)\right)\right]$$ $$=\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-\sum_{i=1}^{N}{x_iE\left(w_i\right)}\right)\left(\sum_{i=1}^{N}{w_iy_i}-\sum_{i=1}^{N}{y_iE\left(w_i\right)}\right)\right]$$ $$=\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-\frac{n}{N}\sum_{i=1}^{N}x_i\right)\left(\sum_{i=1}^{N}{w_iy_i}-\frac{n}{N}\sum_{i=1}^{N}y_i\right)\right]$$ $$=\ E\left[\left(\sum_{i=1}^{N}{w_ix_i}-n\mu_x\right)\left(\sum_{i=1}^{N}{w_iy_i}-n\mu_y\right)\right]$$ $$=\ E\left[\sum_{i=1}^{N}{w_ix_i}\sum_{i=1}^{N}{w_iy_i}-\ n\mu_y\sum_{i=1}^{N}{w_ix_i}-\ n\mu_x\sum_{i=1}^{N}{w_iy_i}+n^2\mu_x\mu_y\right]$$ $$=\ E\left[\sum_{i=1}^{N}{w_ix_iw_iy_i}\right]+E\left[\sum_{i=1}^{N}\sum_{j\neq i}^{N}{w_ix_iw_jy_j}\right]-\ n\mu_y\frac{n}{N}\sum_{i=1}^{N}x_i-\ n\mu_x\frac{n}{N}\sum_{i=1}^{N}y_i+n^2\mu_x\mu_y$$ $$=\ \left(\sum_{i=1}^{N}{x_iy_i}\right)\left(\frac{n}{N}\right)+\left(\sum_{i=1}^{N}\sum_{j\neq i}^{N}{x_iy_j}\right)\left(\frac{n}{N}\right)\left(\frac{n-1}{N-1}\right)-\ n^2\mu_y\mu_x-\ n^2\mu_x\mu_y+n^2\mu_x\mu_y$$ $$=\left(\frac{n\left(N-1\right)}{N\left(N-1\right)}\right)\left(\sum_{i=1}^{N}{x_iy_i}\right)+\left(\frac{n\left(n-1\right)}{N\left(N-1\right)}\right)\left(\sum_{i=1}^{N}\sum_{j\neq i}^{N}{x_iy_j}\right)-n^2\mu_x\mu_y$$ $$=\frac{n\left(N-n\right)}{N\left(N-1\right)}\sum_{i=1}^{N}{x_iy_i}+\frac{n\left(n-1\right)}{N\left(N-1\right)}\sum_{i=1}^{N}{x_iy_i}+\frac{n\left(n-1\right)}{N\left(N-1\right)}\sum_{i=1}^{N}\sum_{j\neq i}^{N}{x_iy_j}-n^2\mu_x\mu_y$$ $$=\frac{n\left(N-n\right)}{\left(N-1\right)}\mu_{xy}+\frac{Nn\left(n-1\right)}{N^2\left(N-1\right)}\sum_{i=1}^{N}x_i\sum_{j=i}^{N}y_j-n^2\mu_x\mu_y$$ $$=\frac{n\left(N-n\right)}{\left(N-1\right)}\mu_{xy}+\frac{Nn\left(n-1\right)}{\left(N-1\right)}\mu_x\mu_y-\frac{n^2\left(N-1\right)}{\left(N-1\right)}\mu_x\mu_y$$ $$=\frac{n\left(N-n\right)}{\left(N-1\right)}\mu_{xy}-\frac{n\left(N-n\right)}{\left(N-1\right)}\mu_x\mu_y$$ $$=\frac{n\left(N-n\right)}{\left(N-1\right)}\left(\mu_{xy}-\mu_x\mu_y\right)$$ $$=\frac{n\left(N-n\right)}{\left(N-1\right)}Cov(x_i,y_i)$$