r'''Calculates the covariance matrix of weighted samples (like the output of an
    importance-sampling run).

    :param samples:

        Matrix-like numpy array; the samples to be used.

    :param weights:

        Vector-like numpy array; the (unnormalized) importance weights.

def calculate_covariance(samples, weights):
    r'''Calculates the covariance matrix of weighted samples (like the output of an
    importance-sampling run).

    :param samples:

        Matrix-like numpy array; the samples to be used.

    :param weights:

        Vector-like numpy array; the (unnormalized) importance weights.

    '''
    assert len(samples) == len(weights), "The number of samples (got %i) must equal the number of weights (got %i)." % (len(samples),len(weights))

    sum_weights_sq = (weights.sum())**2
    sum_sq_weights = (weights**2).sum()

    mean  = calculate_mean(samples, weights)

    return sum_weights_sq / (sum_weights_sq - sum_sq_weights)  *\
           calculate_expectation(samples, weights, lambda x: _np.einsum('i,j', x - mean, x - mean))