# Finding the covariance of two random sums

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.

• Neither proposed answer can possibly be correct, because when $n=N,$ both sums are constant, whence their covariance is zero.
– whuber
Commented Jul 16, 2020 at 20:53
• Thanks for the reply. I didn't intend to suggest that n = N here. But maybe I stated something incorrectly. Where did you infer n=N from? Commented Jul 16, 2020 at 21:05
• I didn't: I offer it as a counterexample. You haven't excluded $n=N$ so presumably you suppose your algebra applies to that case. Since it produces as obviously incorrect answer in that case, your formula cannot hold in general. Any generally valid formula must give $0$ when you plug $N$ in for $n.$
– whuber
Commented Jul 16, 2020 at 21:07
• Good point. Thanks for the clarification. I believe the second equation works with a finite population correction applied. I ran simulations in R that were just about spot on, but that was with a small n relative to N. So I think the last result above basically applies to large N with small n. Does that sound right? (I have edited my post to reflect this.) Commented Jul 16, 2020 at 21:22
• I recommend you redo your simulation with really tiny $N$ so you can detect the difference between $N$ and $N-1$ in the denominator. Note that the correct meaning of "$\operatorname{Cov}$" in this context is the population covariance, not the bias-corrected covariance estimator (there isn't anything to estimate).
– whuber
Commented Jul 16, 2020 at 21:54

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...
$$$$
`
• Thanks for this great answer. There's just one thing I'm having trouble following though. How do you go from: $\ - \frac{n(N-n)}{N^2(N-1)}\sum_{i\ne j}^N x_iy_j \\$ to $\ - \frac{n(N-n)}{N-1} \bar{x}\bar{y} + \frac{n(N-n)}{N(N-1)}\overline{xy}\\$ Commented Jul 17, 2020 at 4:03
• It's just algebra. I think of it as adding and subtracting $\sum_i x_i y_i$ to $\sum_{i\ne j} x_iy_j.$ Adding in this sum gives $\sum_{i,j}x_iy_j=\left(\sum_ix_i\right)\left(\sum_jy_j\right),$ which is a multiple of $\bar{x}\bar{y}.$
– whuber
Commented Jul 17, 2020 at 14:37
• Thanks for the great and insightful answer. I was wondering what if $n$ itself is a binomial random variable, i.e., $n \sim B(N, p)$? Commented Oct 7, 2023 at 10:35
• @Zichao $n$ is fixed: that's why the sum at the beginning has zero variance. You can confirm that by inspecting the code at the end following the comment "Specify the sample size."
– whuber
Commented Oct 7, 2023 at 11:07
• I understand that in this case $n$ is fixed. I just encountered a very similar problem recently where $n$ was a binomial random variable rather than a constant. Could we still derive the covariance of two sums under this condition? Commented Oct 7, 2023 at 11:31

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)$$