python - 使用缺失值计算 scipy 中的成对距离

Question

我对 scipy.spatial.distance.pdist 如何处理缺失的 (nan) 值感到有点困惑。

因此，以防万一我弄乱了矩阵的维度，让我们把它移开。来自文档:

The points are arranged as m n-dimensional row vectors in the matrix X.

那么让我们在 10 维空间中生成三个具有缺失值的点:

numpy.random.seed(123456789)
data = numpy.random.rand(3, 10) * 5
data[data < 1.0] = numpy.nan

如果我计算这三个观测值的欧氏距离:

pdist(data, "euclidean")

我得到:

array([ nan,  nan,  nan])

However, if I filter all the columns with missing values I do get proper distance values:

valid = [i for (i, col) in enumerate(data.T) if ~numpy.isnan(col).any()]
pdist(data[:, valid], "euclidean")

我得到:

array([ 3.35518662,  2.35481185,  3.10323893])

This way, I throw away more data than I'd like since I don't need to filter the whole matrix but only the pairs of vectors being compared at a time. Can I make pdist or a similar function perform pairwise masking, somehow?

---

Edit:

Since my full matrix is rather large, I did some timing tests on the small data set provided here.

1.) The scipy function.

%timeit pdist(data, "euclidean")

10000 loops, best of 3: 24.4 µs per loop

2.) Unfortunately, the solution provided so far is roughly 10 times slower.

%timeit numpy.array([pdist(data[s][:, ~numpy.isnan(data[s]).any(axis=0)], "euclidean") for s in map(list, itertools.combinations(range(data.shape[0]), 2))]).ravel()

1000 loops, best of 3: 231 µs per loop

3.) Then I did a test of "pure" Python and was pleasantly surprised:

from scipy.linalg import norm

%%timeit
m = data.shape[0]
dm = numpy.zeros(m * (m - 1) // 2, dtype=float)
mask = numpy.isfinite(data)
k = 0
for i in range(m - 1):
    for j in range(i + 1, m):
        curr = numpy.logical_and(mask[i], mask[j])
        u = data[i][curr]
        v = data[j][curr]
        dm[k] = norm(u - v)
        k += 1

10000 loops, best of 3: 98.9 µs per loop

所以我认为前进的方向是在函数中对上述代码进行 Cythonize。

Answer 1

如果我理解正确，您需要两个向量具有有效值的所有维度的距离。

不幸的是 pdist 不理解这种意义上的掩码数组，因此我修改了您的半解决方案以不减少信息。然而，它不是最有效的解决方案，也不是最易读的:

np.array([pdist(data[s][:, ~numpy.isnan(data[s]).any(axis=0)], "euclidean") for s in map(list, itertools.combinations(range(data.shape[0]), 2))]).ravel()

把它变成一个数组的外层，ravel 只是为了让它的形状与你期望的相匹配。

itertools.combinations 生成 data 数组的所有成对可能索引。

然后我只是对这些数据进行切片(必须是 list 而不是 tuple 才能正确切片)并对 nan 进行成对过滤> 正如您的代码所做的那样。