用总体平均数计算协方差的NumPy向量化方法（用于调查数据）

Question

我先来介绍一下调查数据。你可以跳过它

简介

调查数据是由调查专家提出的复杂抽样模式形成的。样本可以按州、区或县、地区、地方等对一个国家进行分层。它甚至可以按种族、收入等对人进行分层。一旦创建了地层来解决调查的设计问题，就从这些地层中随机抽取样本。我们对这些样本进行了调查，但由于有些人不喜欢透露他们认为是个人的信息，所以并不是所有的回答都是有效的。为了克服这一点，通常会对地层进行过采样，因此您希望正确有效的响应数量至少等于您所需的样本量。在

当把这些反应放在一起进行分析时，需要对这些反应进行加权，这样分析就可以根据样本来估计整个人群的样子。由于抽样本身非常复杂，而且样本本身是基于其他调查的，因此也存在着调查中的缺点，所以权重不是一个单一的数字。权重有自己的分布。在

如果只想计算数据的中心趋势，比如平均值或中位数，那么最终的权重就足够了。但是，如果没有集中趋势扩散的概念，结果是不完整的。价差是从方差中获得的。对于调查数据，有两个变量在发挥作用，真实的总体方差和由权重引起的方差。这使得方差计算非常复杂，而且容易出错

为了克服方差计算的复杂性，目前有各种各样的渐近方法被广泛应用。其中最受欢迎的是刀切法。Jackknife方法包括删除一部分响应，并重新加权其余响应。因此，一个100阶的折刀将删除1/100的响应，并重新加权剩余的99/100个响应。这将重复100次，每次删除响应的不同1/100部分。完成后，将有一个100刀重（正式称为复制重量）。在

100个重复权重中的每一个，当与总样本权重一样使用时，都将给出一个关于关注的中心趋势的度量。现在，总中心趋势的方差可以从从所有刀切权重得到的中心趋势中计算出来。在

问题陈述

我从总样本权重中计算出总体均值的真实值。我还使用所有的复制权重计算了相同的度量。现在我需要估计方差和协方差，为此，我需要从观测值的平均值中减去观测值。在

在我的例子中，观察值是从重复权重中获得的，平均值是从总权重中获得的。预先打包的协方差函数，如numpy.cov()计算来自传递的数据数组的平均值。在我的情况下，这会给出不正确的答案。基本上，我需要以下伪代码的矢量化版本：

def cov_var(arr, means):
    """
    arr is a matrix of experiments along rows and observations along columns
    means is a vector as long as the number of rows in arr
    """

    for i in range(len(arr)):
        s = 0
        for j in range(len(arr)):
            for k in range(len(arr[i])):
                s += (arr[i][k] - means[i]) * (arr[j][k] - means[j])

Answer 1

这篇文章建议对^{}中列出的矢量化实现进行改进。除了通常的内务错误检查代码，其关键似乎包含在两行-

arr1 = np.subtract(arr.T, means).T
cov = np.dot(arr1, arr1.T)

---

改进1

第一个np.subtract(arr.T, means).T涉及两个转置和一个函数调用。这可以通过用^{}扩展means的维度并执行减法来实现。现在，这两种方法都将使用broadcasting，这基本上就是在这里完成的计算，因此我不期望它有任何显著的加速。替代守则将是-

^{pr2}$

运行时测试和验证-

In [16]: arr = np.random.rand(5000,5000)
    ...: means = np.random.rand(5000)
    ...: 

In [17]: np.allclose(arr - means[:,None],np.subtract(arr.T, means).T)
Out[17]: True

In [18]: %timeit np.subtract(arr.T, means).T
1 loops, best of 3: 346 ms per loop

In [19]: %timeit arr - means[:,None]
1 loops, best of 3: 332 ms per loop

改进2

接下来，np.dot(arr1, arr1.T)可以用np.einsum实现，如下-

np.einsum('ij,kj',arr1,arr1)

运行时测试和验证-

In [16]: arr = np.random.rand(200,300)
    ...: means = np.random.rand(200)
    ...: 

In [17]: arr1 = np.subtract(arr.T, means).T

In [18]: np.allclose(np.dot(arr1, arr1.T),np.einsum('ij,kj',arr1,arr1))
Out[18]: True

In [19]: %timeit np.dot(arr1, arr1.T)
100 loops, best of 3: 17.8 ms per loop

In [20]: %timeit np.einsum('ij,kj',arr1,arr1)
100 loops, best of 3: 9.12 ms per loop

Answer 2

经过一番思考和几次测试，我想出了以下简单的代码：

def var_covar(arr, means=None, axis=1, normalize=True, bias=None):
    """
    Computes the covariance of array arr. Optionally it takes means
    separately instead of computing it from the data.
         
    PARAMETERS
    :arr: array_like; A NumPy array or other data type that can be converted
        into an NumPy array.
    :means: Optional, array_like; A vector of means for each observation in
        arr. Mostly would be population means.
    :axis: Optional, int, default 1; Operates either by row or by column.
        Default assumes experiments are along rows, and observations are
        along columns.
    :normalize: Optional, bool, default True; Normalizes the data by 1/(N-bias),
        where N is number of observations.
    :bias: Optional, float; Bias is subtracted from N, to account for sample
        mean. If means are passed, it is assumed they are population means
        and so bias is set to 0. Else to 1. Passing a value overrides the
        default value of bias.
         
    RETURNS
    :cov: array_likep; A 2-D NumPy array if arr is a matrix. A float if arr
        is a vector
         
    NOTES:
      Maximum number of supported dimensions is 2
      It is assumed that passed means are population means, and not sample
      means. This assumption affects bias. Pass a value for bias to override
      the default bias.
      It is recommended that length of means equals number of observations in
      arr. If that is not the case, means will be broadcast, but it may throw
      an error.
      axis must be either 0 or 1. Higher dimensions are not supported.
    """

    arr = np.asarray(arr)
    tbias = bias

    if arr.ndim > 2:
        raise RuntimeError('Cannot handle arr larger than a 2x2 matrix')

    if axis not in [0,1]:
        raise ValueError('axis must be either 0 or 1')
    elif axis == 0:
        arr = arr.T

    if means is None:
        means = np.average(arr, axis=1)
        bias = 1
    else:
        bias = 0
    means = np.asarray(means)

    arr = np.subtract(arr.T, means).T
    cov = np.dot(arr, arr.T)

    if normalize:
        if tbias is not None:
            bias = tbias
        factor = 1/(arr.shape[0] - bias)
        cov = np.multiply(cov, factor)

    return cov

它比较快，完全矢量化。我把这个分享给其他需要它的人。在

如果你认为你做得比我好，请随时发表答案。我很乐意接受它作为答案。在