np.random.seed(987467)
x = np.random.multivariate_normal([0, 1.0], [[1.0, 0.5], [0.5, 1]], size=20)
weights = np.random.randint(1, 4, size=20)
xlong = np.repeat(x, weights, axis=0)
ds = DescrStatsW(x, weights=weights)
print('cov statsmodels')
print(ds.cov)
self = ds
ds_cov = np.dot(self.weights * self.demeaned.T, self.demeaned
    ) / self.sum_weights
print('\nddof=0')
print(ds_cov)
print(np.cov(xlong.T, bias=1))
ds_cov0 = np.dot(self.weights * self.demeaned.T, self.demeaned) / (self.
    sum_weights - 1)
print('\nddof=1')
print(ds_cov0)
print(np.cov(xlong.T, bias=0))
Python package that supports weighted covariance computation