Correct equation for weighted unbiased sample covariance

Question

I'm looking for the correct equation to compute the weighted unbiased sample covariance. Internet sources are quite rare on this theme and they all use different equations.

The most likely equation I've found is this one:

$q_{jk}=\frac{\sum_{i=1}^{N}w_i}{\left(\sum_{i=1}^{N}w_i\right)^2-\sum_{i=1}^{N}w_i^2} \sum_{i=1}^N w_i \left( x_{ij}-\bar{x}_j \right) \left( x_{ik}-\bar{x}_k \right) .$

From: https://en.wikipedia.org/wiki/Sample_mean_and_sample_covariance#Weighted_samples

Of course, you have to compute the weighted (unbiased) sample mean beforehand.

However, I have found several other formulas like:

$q_{jk}= \frac{1}{\sum_{i=1}^N w_i)-1}\sum_{i=1}^N w_i \left( x_{ij}-\bar{x}_j \right) \left( x_{ik}-\bar{x}_k \right) .$

Or I've even seen some source codes and academic papers just using the standard covariance formula but with the weighted sample mean instead of the sample mean...

Can someone help me and shed some light?

/EDIT: my weights are simply the number of observations for a sample in the dataset, thus weights.sum() = n

Answer 1

Found the solution in a 1972's book (George R. Price, Ann. Hum. Genet., Lond, pp485-490, Extension of covariance selection mathematics, 1972).

Biased weighted sample covariance:

$\Sigma=\frac{1}{\sum_{i=1}^{N}w_i}\sum_{i=1}^N w_i \left(x_i - \mu^*\right)^T\left(x_i - \mu^*\right)$

And the unbiased weighted sample covariance given by applying the Bessel correction:

$\Sigma=\frac{1}{\sum_{i=1}^{N}w_i - 1}\sum_{i=1}^N w_i \left(x_i - \mu^*\right)^T\left(x_i - \mu^*\right)$

Where $\mu^*$ is the (unbiased) weighted sample mean:

$\mathbf{\mu^*}=\frac{\sum_{i=1}^N w_i \mathbf{x}_i}{\sum_{i=1}^N w_i}$

Important Note: this works only if the weights are "repeat"-type weights, meaning that each weight represent the number of occurrences of one observation, and that $\sum_{i=1}^N w_i=N^*$ where $N^*$ represent the real sample size (real total number of samples, accounting for the weights).

I have updated the article on Wikipedia, where you will also find the equation for unbiased weighted sample variance:

https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Weighted_sample_covariance

Practical note: I advise you to first multiply column-by-column $w_i$ and $\left(x_i - \mu^*\right)$ and then do a matrix multiplication with $\left(x_i - \mu^*\right)$ to wrap things up and automatically perform the summation. Eg in Python Pandas/Numpy code:

import pandas as pd
import numpy as np
# X is the dataset, as a Pandas' DataFrame
mean = mean = np.ma.average(X, axis=0, weights=weights) # Computing the weighted sample mean (fast, efficient and precise)
mean = pd.Series(mean, index=list(X.keys())) # Convert to a Pandas' Series (it's just aesthetic and more ergonomic, no differenc in computed values)
xm = X-mean # xm = X diff to mean
xm = xm.fillna(0) # fill NaN with 0 (because anyway a variance of 0 is just void, but at least it keeps the other covariance's values computed correctly))
sigma2 = 1./(weights.sum()-1) * xm.mul(weights, axis=0).T.dot(xm); # Compute the unbiased weighted sample covariance

Did a few sanity checks using a non-weighted dataset and an equivalent weighted dataset, and it works correctly.

For more details about the theory of unbiased variance/covariance, see this post.

Answer 2

Actually, pandas is not required here, outer product from numpy will suffice and should be faster to compute the weighted (unbiased) sample covariance matrix, using just couple of lines of code:

mu = np.sum(weights.reshape(-1,1)*X, axis=0) / np.sum(weights)
Sigma2 = (weights.reshape(-1,1)*(X-mu)).T@(X-mu) / (weights.sum()-1)

Here $X$ is the data matrix (numpy ndarray) and weights is the vector of weights, one for each observation (row) in $X$.

We can test from the next code snippet that it produces the same result and is faster (at least with a small synthetic data):

with numpy only:

from time import time
import numpy as np
assuming that the same data is stored as matrix
np.random.seed(1)
X1 = np.random.randint(0,100,size=(100, 4))
weights = np.random.randint(0,100,size=100)
start = time()
mu = np.sum(weights.reshape(-1,1)X1, axis=0) / np.sum(weights)
Sigma2 = (weights.reshape(-1,1)(X1-mu)).T@(X1-mu) / (weights.sum()-1)
print('time taken', time() - start)
Sigma2
time taken 0.00099945068359375
array([[ 890.94399402,  -85.64453306,  115.59965823,   27.03209622],
[ -85.64453306,  930.63865593, -219.27393837, -147.81191515],
[ 115.59965823, -219.27393837,  752.47415657,   92.31686826],
[  27.03209622, -147.81191515,   92.31686826,  745.7993326 ]])

with pandas and numpy:

import pandas as pd
# assuming that the data is stored as dataframe
X = pd.DataFrame(X1, columns=list('ABCD'))
weights = np.random.randint(0,100,size=100)
start = time()
mean = np.ma.average(X, axis=0, weights=weights) # Computing the weighted sample mean (fast, efficient and precise)
mean = pd.Series(mean, index=list(X.keys())) # Convert to a Pandas' Series (it's just aesthetic and more ergonomic, no differenc in computed values)
xm = X-mean # xm = X diff to mean
xm = xm.fillna(0) 
sigma2 = 1./(weights.sum()-1) * xm.mul(weights, axis=0).T.dot(xm); 
print('time taken', time() - start)
sigma2
time taken 0.0059969425201416016
A B   C   D
A 890.943994  -85.644533  115.599658  27.032096
B -85.644533  930.638656  -219.273938 -147.811915
C 115.599658  -219.273938 752.474157  92.316868
D 27.032096   -147.811915 92.316868   745.799333