我正在寻找计算两个 numpy 数组(x 和 y)之间马氏距离的 NumPy 方法。 以下代码可以使用 Scipy 的 cdist 函数正确计算相同的内容。由于这个函数在我的例子中计算了不必要的矩阵,我想要更直接的方法来仅使用 NumPy 来计算它。
import numpy as np
from scipy.spatial.distance import cdist
x = np.array([[[1,2,3,4,5],
[5,6,7,8,5],
[5,6,7,8,5]],
[[11,22,23,24,5],
[25,26,27,28,5],
[5,6,7,8,5]]])
i,j,k = x.shape
xx = x.reshape(i,j*k).T
y = np.array([[[31,32,33,34,5],
[35,36,37,38,5],
[5,6,7,8,5]],
[[41,42,43,44,5],
[45,46,47,48,5],
[5,6,7,8,5]]])
yy = y.reshape(i,j*k).T
results = cdist(xx,yy,'mahalanobis')
results = np.diag(results)
print results
[ 2.28765854 2.75165028 2.75165028 2.75165028 0. 2.75165028
2.75165028 2.75165028 2.75165028 0. 0. 0. 0.
0. 0. ]
我的试用:
VI = np.linalg.inv(np.cov(xx,yy))
print np.sqrt(np.dot(np.dot((xx-yy),VI),(xx-yy).T))
谁能纠正这个方法?
公式如下:
最佳答案
我认为你的问题在于你的协方差矩阵的构造。尝试:
X = np.vstack([xx,yy])
V = np.cov(X.T)
VI = np.linalg.inv(V)
print np.diag(np.sqrt(np.dot(np.dot((xx-yy),VI),(xx-yy).T)))
输出:
[ 2.28765854 2.75165028 2.75165028 2.75165028 0. 2.75165028
2.75165028 2.75165028 2.75165028 0. 0. 0. 0.
0. 0. ]
要在此处隐式创建中间数组的情况下执行此操作,您可能不得不为 Python 循环牺牲 C 循环:
A = np.dot((xx-yy),VI)
B = (xx-yy).T
n = A.shape[0]
D = np.empty(n)
for i in range(n):
D[i] = np.sqrt(np.sum(A[i] * B[:,i]))
编辑:实际上,使用 np.einsum
voodoo,您可以删除 Python 循环并大大加快它的速度(在我的系统上,从 84.3 µs 到 2.9 µs):
D = np.sqrt(np.einsum('ij,ji->i', A, B))
编辑:正如@Warren Weckesser 指出的那样,einsum
也可用于消除中间A
和B
数组:
delta = xx - yy
D = np.sqrt(np.einsum('nj,jk,nk->n', delta, VI, delta))
关于python - 仅使用 NumPy 计算马氏距离,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/27686240/