用 Mahalanobis 处理奇异矩阵答案

【问题标题】：Handling singular matrix with Mahalanobis用 Mahalanobis 处理奇异矩阵
【发布时间】：2021-01-20 16:13:47
【问题描述】：

我有一个组数据框，我想计算每个组的马氏距离。

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

我的数据如下：

foo <- data.frame(GRP = c("a","a","b","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 的建议，我在auto.mahalanobis 函数中添加了tol 参数，以避免在计算小数的协方差矩阵时出现问题。

然后我尝试将此函数与我的数据集一起使用，并收到以下关于奇异矩阵的错误：

 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".

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

【问题讨论】：

标签： r dplyr mahalanobis

【解决方案1】：

最简单的方法可能是：

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)；我认为这可能（？）可以作为替代品，但尚未检查数学。

【讨论】：