r - R中的马氏距离

Question

我在包 StatMatch ( http://cran.r-project.org/web/packages/StatMatch/StatMatch.pdf ) 中找到了 mahalanobis.dist 函数，但它并没有完全符合我的要求。它似乎正在计算从 data.y 中的每个观察到 data.x 中的每个观察的马氏距离

我想计算 data.y 中的一个观测值与 data.x 中的所有观测值的马氏距离。如果有意义的话，基本上计算一个点到点“云”的马氏距离。有点理解一个观察是另一组观察的一部分的概率

这个人（http://people.revoledu.com/kardi/tutorial/Similarity/MahalanobisDistance.html）似乎正在这样做，我试图在R中复制他的过程，但是当我到达底部时它失败了等式的：

mahaldist = sqrt((inversepooledcov %*% t(meandiffmatrix)) %*% meandiffmatrix)

我正在使用的所有代码都在这里：

a = rbind(c(2,2), c(2,5), c(6,5),c(7,3))

colnames(a) = c('x', 'y')

b = rbind(c(6,5),c(3,4))

colnames(b) = c('x', 'y')

acov = cov(a)
bcov = cov(b)

meandiff1 = mean(a[,1]) - mean(b[,1])

meandiff2 = mean(a[,2]) - mean(b[,2])

meandiffmatrix = rbind(c(meandiff1,meandiff2))

totaldata = dim(a)[1] + dim(b)[1]

pooledcov = (dim(a)[1]/totaldata * acov) + (dim(b)[1]/totaldata * bcov)

inversepooledcov = solve(pooledcov)

mahaldist = sqrt((inversepooledcov %*% t(meandiffmatrix)) %*% meandiffmatrix)

Answer 1

如何使用包mahalanobis中的功能stats：

 mahalanobis(x, center, cov, inverted = FALSE, ...)

Answer 2

我一直在您查看的同一个网站上尝试过这个，然后偶然发现了这个问题。我设法让脚本工作，但我得到了不同的结果。

#WORKING EXAMPLE
#MAHALANOBIS DIST OF TWO MATRICES

#define matrix
mat1<-matrix(data=c(2,2,6,7,4,6,5,4,2,1,2,5,5,3,7,4,3,6,5,3),nrow=10)
mat2<-matrix(data=c(6,7,8,5,5,5,4,7,6,4),nrow=5)
#center data
mat1.1<-scale(mat1,center=T,scale=F)
mat2.1<-scale(mat2,center=T,scale=F)
#cov matrix
mat1.2<-cov(mat1.1,method="pearson")
mat2.2<-cov(mat2.1,method="pearson")
n1<-nrow(mat1)
n2<-nrow(mat2)
n3<-n1+n2
#pooled matrix
mat3<-((n1/n3)*mat1.2) + ((n2/n3)*mat2.2)
#inverse pooled matrix
mat4<-solve(mat3)
#mean diff
mat5<-as.matrix((colMeans(mat1)-colMeans(mat2)))
#multiply
mat6<-t(mat5) %*% mat4
#multiply
sqrt(mat6 %*% mat5)

我认为该函数mahalanobis()用于计算一个矩阵中个体（行）之间的马氏距离。pairwise.mahalanobis()from函数package(HDMD)可以比较两个或多个矩阵并给出矩阵之间的马氏距离。

Answer 3

您可以将函数包装stats::mahalanobis如下以输出马氏距离矩阵（成对马氏距离）：

# x - data frame
# cx - covariance matrix; if not provided, 
#      it will be estimated from the data
mah <- function(x, cx = NULL) {
  if(is.null(cx)) cx <- cov(x)
  out <- lapply(1:nrow(x), function(i) {
    mahalanobis(x = x, 
                center = do.call("c", x[i, ]),
                cov = cx)
  })
  return(as.dist(do.call("rbind", out)))
}

然后，您可以对数据进行聚类并绘制它，例如：

# Dummy data
x <- data.frame(X = c(rnorm(10, 0), rnorm(10, 5)), 
                Y = c(rnorm(10, 0), rnorm(10, 7)), 
                Z = c(rnorm(10, 0), rnorm(10, 12)))
rownames(x) <- LETTERS[1:20]
plot(x, pch = LETTERS[1:20])

# Comute the mahalanobis distance matrix
d <- mah(x)
d

# Cluster and plot
hc <- hclust(d)
plot(hc)

Answer 4

取平方根之前的输出是：

inversepooledcov %*% t(meandiffmatrix) %*% meandiffmatrix
          [,1]        [,2]
x -0.004349227 -0.01304768
y  0.114529639  0.34358892

我认为你不能取负数的平方根，所以你有NAN负元素：

 sqrt(inversepooledcov %*% t(meandiffmatrix) %*% meandiffmatrix)
       [,1]      [,2]
x       NaN       NaN
y 0.3384223 0.5861646

Warning message:
In sqrt(inversepooledcov %*% t(meandiffmatrix) %*% meandiffmatrix) :
  NaNs produced

Answer 5

如果协方差矩阵是恒等式，则马氏距离等于（平方）欧几里得距离。如果变量之间存在协方差，则可以通过首先对矩阵进行白化以消除协方差来使 Mahalanobis 和 sq Euclidean 相等。即，做：

#X is your matrix
if (!require("whitening")) install.packages("whitening")

X <- whitening::whiten(X) # default is ZCA (Mahalanobis) whitening
X_dist <- dist(X, diag = T, method = "euclidean")^2

您可以确认这为您提供了与 Davit 在先前答案之一中提供的代码相同的距离矩阵。

Answer 6

使用 R 包“biotools”有一个非常简单的方法。在这种情况下，您将得到一个平方距离马氏矩阵。

#Manly (2004, p.65-66)

x1 <- c(131.37, 132.37, 134.47, 135.50, 136.17)
x2 <- c(133.60, 132.70, 133.80, 132.30, 130.33)
x3 <- c(99.17, 99.07, 96.03, 94.53, 93.50)
x4 <- c(50.53, 50.23, 50.57, 51.97, 51.37)

#size (n x p) #Means 
x <- cbind(x1, x2, x3, x4) 

#size (p x p) #Variances and Covariances
Cov <- matrix(c(21.112,0.038,0.078,2.01, 0.038,23.486,5.2,2.844, 
        0.078,5.2,24.18,1.134, 2.01,2.844,1.134,10.154), 4, 4)

library(biotools)
Mahalanobis_Distance<-D2.dist(x, Cov)
print(Mahalanobis_Distance)

Answer 7

您现在可以通过metan包计算马氏距离。参考函数mahala()和mahala_design(). 包装文件