python - 用于计算总体平均值协方差的 NumPy 矢量化方法(针对调查数据)

Question

我将首先介绍一下调​​查数据。您可以跳过它

简介

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

将响应放在一起进行分析时，需要对响应进行加权，以便分析根据样本估计整个总体的外观。由于抽样本身非常复杂，并且样本本身是基于其他调查的，因此受到该调查的缺点的影响，因此权重不是一个单一的数字。权重有自己的分布。

如果只想计算数据的集中趋势，例如平均值或中位数，则最终权重应该足够了。但如果不考虑集中趋势的传播，结果是不完整的。分布是从方差获得的。对于调查数据，存在两种方差:真实总体方差和权重引起的方差。这使得方差计算极其复杂且容易出错

当今有多种渐近方法被广泛使用来克服方差计算的复杂性。最流行的方法之一是折刀法。 Jackknife 方法涉及删除一部分响应，并对剩余响应重新加权。因此，100 阶 Jackknife 将删除 1/100 的响应，并重新权衡剩余的 99/100 的响应。这将重复 100 次，每次删除响应的不同 1/100 部分。完成后，您将拥有 100 个折刀配重(正式称为复制配重)。

当像总样本权重一样使用时，100 个重复权重中的每一个都会给出感兴趣的集中趋势的度量。现在，可以根据所有 Jackknife 权重获得的集中趋势来计算总集中趋势的方差。

问题陈述

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

在我的例子中，观察值是从重复权重中获得的，平均值是从总重量中获得的。预先打包的协方差函数，如 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

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

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

---

改进#1

第一个np.subtract(arr.T,means).T涉及两次转置和一次函数调用。这可以通过使用 None/np.newaxis 扩展 means 的维度来简单地实现。并进行减法。现在，这两种方法都将使用广播，这本质上就是这里完成的计算，因此我不期望它有任何显着的加速。替代代码是 -

arr1 = arr - means[:,None]

运行时测试和验证 -

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