用Mahalanobis处理奇异矩阵

Question

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

我正在尝试将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".

同样的功能适用于样本大小更大的其他组，有没有建议的方法来解决这个问题，或者在样本太小时跳过这些组？

Answer 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)))。

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