用 Mahalanobis 处理奇异矩阵

问题描述

我有一个组数据框，我想计算每个组的马哈拉诺比斯距离。

我正在尝试将 Mahalanobis 函数应用于数百个组，但由于样本量小（只有两行），一个特定的组导致了问题。

我的数据如下所示：

foo <- data.frame(GRP = c("a","a","b","b"),X = c(1,1,15,12,50),Y = c(2.17,12.44,50,70,100))

我从here那里借用了一个函数思想，它看起来如下：

auto.mahalanobis <- function(temp) {
 mahalanobis(temp,center = colMeans(temp,na.rm=T),cov = cov(temp,use="pairwise.complete.obs"),tol=1e-20,inverted = FALSE
             )
}

根据建议 here，我向 tol 函数添加了 auto.mahalanobis 参数，以避免在计算小数协方差矩阵时出现问题。

然后我尝试在我的数据集上使用这个函数，但出现以下关于奇异矩阵的错误：

 z <- foo %>% group_by(GRP) %>% mutate(mahal = auto.mahalanobis(data.frame(X,Y)))

Error: Problem with `mutate()` input `mahal`.
x Lapack routine dgesv: system is exactly singular: U[1,1] = 0
i Input `mahal` is `auto.mahalanobis(data.frame(X,Y))`.
i The error occurred in group 1: GRP = "a".

相同的功能适用于样本量较大的其他组，是否有建议的方法来解决此问题或在样本太小时跳过此类组？

解决方法

最简单的方法可能是：

auto.mahalanobis <- function(temp) {
 m <- try(silent=TRUE,mahalanobis(temp,center = colMeans(temp,na.rm=TRUE),cov = cov(temp,use="pairwise.complete.obs"),tol=1e-20,inverted = FALSE
             ))
 if (!inherits(m,"try-error")) return(m)
 return(rep(NA_real_,length(temp))
}

（未经测试：真正的程序员可能会改用 tryCatch()）

如果您认为问题只会在 n==2 时发生，您可以使用 if 子句，例如if (length(temp)<min_length) return(rep(NA_real_,length(temp)))。

或者，您可以制作一个使用广义逆 (mahalanobis()) 而不是常规矩阵逆 (MASS::ginv) 的黑客版本的 solve；我认为这可能 (?) 可以作为替代品，但还没有检查数学。

dplyr dplyr mahalanobis r r