Correct equation for weighted unbiased sample covariance 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
 A: 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.
