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

我对scipy.spatial.distance.pdist处理missing(nan)值的方式感到有点困惑.

所以,万一我弄乱了矩阵的维度,让我们把它弄清楚.来自文档:

这些点被排列为矩阵X中的m个n维行向量.

因此,让我们在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])

但是,如果我过滤掉所有缺少值的列,我会得到适当的距离值:

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])

这样,我丢弃了比我想要的更多的数据,因为我不需要过滤整个矩阵,而只需要一次比较一对矢量.我可以制作pdist或类似的功能执行成对屏蔽,不知何故？

编辑:

由于我的完整矩阵相当大,我对这里提供的小数据集进行了一些时序测试.

1.)scipy功能.

%timeit pdist(data, "euclidean")

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

2.)不幸的是,到目前为止提供的解决方案大约慢了10倍.

%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.)然后我做了一个"纯粹的"Python测试,并感到惊喜:

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

所以我认为前进的方法是在函数中Cython化上面的代码.