计算经验/样本协方差函数的快速优雅方法

Question

有人知道在Python中计算经验/样本协方差图的好方法吗？

这是一个书籍的截图，其中包含了协方差图的良好定义：

如果我理解正确的话，对于给定的滞后/宽度h，我应该得到所有距离小于等于h的点对，并将其值相乘。对于每个这样的点，计算其平均值，即在这种情况下定义为m(x_i)。然而，根据m(x_{i})的定义，如果我想要计算m(x1)，我需要获得距离x1不超过h的值的平均值。这看起来像是一个非常繁重的计算。

首先，我是否理解正确？如果是这样，假设在二维空间中，有什么好的方法来计算这个问题呢？我试图用Python编写代码（使用numpy和pandas），但需要几秒钟时间，而且我甚至不确定它是否正确，这就是为什么我会避免在此处发布代码的原因。以下是另一种非常天真的实现尝试：

from scipy.spatial.distance import pdist, squareform
distances = squareform(pdist(np.array(coordinates))) # coordinates is a nx2 array
z = np.array(z) # z are the values
cutoff = np.max(distances)/3.0 # somewhat arbitrary cutoff
width = cutoff/15.0
widths = np.arange(0, cutoff + width, width)
Z = []
Cov = []

for w in np.arange(len(widths)-1): # for each width
    # for each pairwise distance
    for i in np.arange(distances.shape[0]): 
        for j in np.arange(distances.shape[1]): 
            if distances[i, j] <= widths[w+1] and distances[i, j] > widths[w]:
                m1 = []
                m2 = []
                # when a distance is within a given width, calculate the means of
                # the points involved
                for x in np.arange(distances.shape[1]):
                    if distances[i,x] <= widths[w+1] and distances[i, x] > widths[w]:
                        m1.append(z[x])

                for y in np.arange(distances.shape[1]):
                    if distances[j,y] <= widths[w+1] and distances[j, y] > widths[w]:
                        m2.append(z[y])

                mean_m1 = np.array(m1).mean() 
                mean_m2 = np.array(m2).mean()
                Z.append(z[i]*z[j] - mean_m1*mean_m2)
    Z_mean = np.array(Z).mean() # calculate covariogram for width w
    Cov.append(Z_mean) # collect covariances for all widths

然而，现在我已经确认我的代码存在错误。我知道这是因为我使用变程图来计算协方差函数（协方差函数(h) = 协方差函数(0) - 变程图(h)），并且得到了不同的图形：

并且它应该看起来像这样：

最后，如果您知道一个Python/R/MATLAB库来计算经验协变函数，请告诉我。至少这样我可以验证我所做的。

- r_31415

---

你所写的 m 方程没有意义。如果你对 i 求和，那么在求和之外用 i 作为索引就没有任何意义（例如，在 m(x_i) 中）；也就是说，右侧没有 i。 - tom10

我认为这并不像你想象的那么难，但如果没有正确的方程式，很难知道该怎么做。我也不想读一整章的书来找出正确的方程式。我会等一段时间，但最终我会投票关闭，因为问题不明确。基本上，你只需要对成对的度量进行求和，这在numpy中可能很容易实现。但首先你需要知道你需要做什么。 - tom10

好的，我可以尝试更改方程式。顺便说一下，如果您点击我之前评论中的链接，它会引导到方程式。请告诉我如何改进问题，而不是投票关闭。 - r_31415

在这个例子中，我认为有150个点。然而，它应该用于成千上万的点。我不确定是否可能在不考虑每对距离的情况下计算协方差函数/变异函数。 - r_31415

太好了！感谢您的评论。 - r_31415

Answer 1

您可以使用scipy.cov，但如果直接进行计算（这很容易），则有更多方法可以加快速度。

首先，通过先制作空间相关性，然后使用根据底图生成的随机数据点来生成数据，其中数据根据底图定位，并且也采用底图的值。

编辑1：
我更改了数据点生成器，因此位置是纯随机的，但z值与空间地图成比例。并且，我更改了地图，以便左侧和右侧相对移动，从而在较大的h处创建负相关。

from numpy import *
import random
import matplotlib.pyplot as plt

S = 1000
N = 900
# first, make some fake data, with correlations on two spatial scales
#     density map
x = linspace(0, 2*pi, S)
sx = sin(3*x)*sin(10*x)
density = .8* abs(outer(sx, sx))
density[:,:S//2] += .2
#     make a point cloud motivated by this density
random.seed(10)  # so this can be repeated
points = []
while len(points)<N:
    v, ix, iy = random.random(), random.randint(0,S-1), random.randint(0,S-1)
    if True: #v<density[ix,iy]:
        points.append([ix, iy, density[ix,iy]])
locations = array(points).transpose()
print locations.shape
plt.imshow(density, alpha=.3, origin='lower')
plt.plot(locations[1,:], locations[0,:], '.k')
plt.xlim((0,S))
plt.ylim((0,S))
plt.show()
#     build these into the main data: all pairs into distances and z0 z1 values
L = locations
m = array([[math.sqrt((L[0,i]-L[0,j])**2+(L[1,i]-L[1,j])**2), L[2,i], L[2,j]] 
                         for i in range(N) for j in range(N) if i>j])

这将产生：

上面只是模拟数据，我没有试图优化它的生成等。我假设这就是OP要开始的地方，因为数据已经存在于实际情况中。

---

现在计算“协方差函数”（比生成假数据容易多了）。这里的想法是按h对所有配对和相关值进行排序，然后使用ihvals对它们进行索引。也就是说，到索引ihval为止的求和是方程中N(h)的总和，因为这包括所有h小于所需值的配对。

编辑2：
如下面的评论所建议的，N(h)现在仅限于h-dh和h之间的配对，而不是0和h之间的所有配对（其中dh是ihvals中h值的间距——即，下面使用了S/1000）。

# now do the real calculations for the covariogram
#    sort by h and give clear names
i = argsort(m[:,0])  # h sorting
h = m[i,0]
zh = m[i,1]
zsh = m[i,2]
zz = zh*zsh

hvals = linspace(0,S,1000)  # the values of h to use (S should be in the units of distance, here I just used ints)
ihvals = searchsorted(h, hvals)
result = []
for i, ihval in enumerate(ihvals[1:]):
    start, stop = ihvals[i-1], ihval
    N = stop-start
    if N>0:
        mnh = sum(zh[start:stop])/N
        mph = sum(zsh[start:stop])/N
        szz = sum(zz[start:stop])/N
        C = szz-mnh*mph
        result.append([h[ihval], C])
result = array(result)
plt.plot(result[:,0], result[:,1])
plt.grid()
plt.show()

对于我来说，这看起来很合理，因为可以看到在预期的 h 值处有凸起或凹槽，但我没有进行仔细的检查。

与 scipy.cov 相比，主要的加速是可以预先计算所有的乘积 zz。否则，每次出现新的 h，都需要将 zh 和 zsh 输入到 cov 中，并重新计算所有乘积。可以通过做部分求和来进一步加快计算速度，即从时间步骤 n 的 ihvals[n-1] 到 ihvals[n]，但我觉得这可能没有必要。

- tom10

---

谢谢！我马上会看一下这个。 - r_31415

谢谢你的回答。你似乎将m_+h和m_-h解释为距离h以下的平均值。也许你是对的，但我将其解释为Cov(h)定义中每个值在距离h内的平均值。显然，我的解释需要更加密集的计算。你知道哪里可以确认你的解释吗？ - r_31415

另一个问题是N(h)的定义。我认为正确的定义是距离h内的一组成对数据，但不包括之前的距离。否则，当我们增加滞后时，协方差函数几乎总是为零，因为此时N(h)包括几乎每一对点。此外，我看到过几个实现中都包括一个阈值，以便捕获在h-e和h+e之间的点，其中e是容差值。 - r_31415

@Robert：我认为你关于 N(h) 的看法是正确的，根据定义，我的累加方式是不正确的，所以总和应该是 zh[ihval-1:ihval] ，其中“容差值”是 h 的跨度。(这个问题是否也回答了 m 和“密集计算”的差异，或者你认为还有其他的不同之处？) 我会尽快修改我的答案。 - tom10

@RobertSmith：我已经进行了编辑，包括一个分布，应该在大h时给出负相关性。但是我仍然不明白你对m的其他解释是什么。我认为我的解释是合理的，因为它类似于平均值通常从方差计算中得出的方式。（另外，如果您有完整的推导，您应该能够检查正确的解释，这应该基于数学而不是周围的文字，或者您可以在统计堆栈交换上询问。） - tom10