python - 快速、优雅的方法来计算经验/样本协方差图

Question

如果可能的话，有人知道用 Python 计算经验/样本协方差图的好方法吗？

这是 book 的截图其中包含协变函数的良好定义:

如果我理解正确，对于给定的滞后/宽度 h，我应该得到所有由 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

但是，现在我已经确认我的代码有错误。我知道这是因为我使用变差函数来计算协变函数 (covariogram(h) = covariogram(0) - variogram(h)) 并且我得到了一个不同的图:

它应该是这样的:

最后，如果您知道用于计算经验协方差图的 Python/R/MATLAB 库，请告诉我。至少，这样我可以验证我做了什么。

Answer 1

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

首先，制作一些具有某些空间相关性的假数据。为此，我将首先建立空间相关性，然后使用由此生成的随机数据点，其中数据根据底层 map 定位，并采用底层 map 的值。

编辑 1:
我更改了数据点生成器，因此位置完全是随机的，但 z 值与空间 map 成正比。而且，我更改了 map ，使左侧和右侧相对于彼此移动，以在大 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 是 h 的间距 - ihvals -- 即下面使用了 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，并且所有的乘积都将被重新计算.通过进行部分求和，即从 ihvals[n-1] 到 ihvals[n] 在每个时间步 n，但我怀疑是否有必要。