【问题标题】:Mahalanobis distance of each pair of observations每对观测值的马氏距离
【发布时间】:2017-04-22 21:09:05
【问题描述】:

我正在尝试计算数据集dat 的每个观测值之间的马氏距离,其中每一行是一个观测值,每一列是一个变量。这样的距离定义为:

我写了一个函数来做,但我觉得它很慢。有没有更好的方法在 R 中计算这个?

生成一些数据来测试功能:

generateData <- function(nObs, nVar){
  library(MASS)
  mvrnorm(n=nObs, rep(0,nVar), diag(nVar))
  }

这是我目前写的函数。它们都有效,对于我的数据(800 个 obs 和 90 个变量),method = "forLoop"method = "apply" 分别需要大约 30 和 33 秒。

mhbd_calc2 <- function(dat, method) { #Method is either "forLoop" or "apply"
  dat <- as.matrix(na.omit(dat))
  nObs <- nrow(dat)
  mhbd <- matrix(nrow=nObs,ncol = nObs)
  cv_mat_inv = solve(var(dat))

  distMH = function(x){  #Mahalanobis distance function
    diff = dat[x[1],]-dat[x[2],]
    diff %*% cv_mat_inv %*% diff
  }

  if(method=="forLoop")
  {
    for (i in 1:nObs){
      for(j in 1:i){
        mhbd[i,j] <- distMH(c(i,j))
      }
    }
  }
  if(method=="apply")
  {
    mhbd[lower.tri(mhbd)] = apply(combn(nrow(dat),2),2, distMH)
  }
  result = sqrt(mhbd)
  colnames(result)=rownames(dat)
  rownames(result)=rownames(dat)
  return(as.dist(result))
}

注意:我尝试使用 outer(),但速度更慢(60 秒)

【问题讨论】:

    标签: r matrix distance mahalanobis


    【解决方案1】:

    你需要一些数学知识。

    1. 对经验协方差进行 Cholesky 分解,然后标准化您的观察结果;
    2. 使用dist 计算变换观测值的欧几里得距离。

    dist.maha <- function (dat) {
      X <- as.matrix(na.omit(dat))  ## ensure a valid matrix
      V <- cov(X)  ## empirical covariance; positive definite
      L <- t(chol(V))  ## lower triangular factor
      stdX <- t(forwardsolve(L, t(X)))  ## standardization
      dist(stdX)  ## use `dist`
      }
    

    示例

    set.seed(0)
    x <- matrix(rnorm(6 * 3), 6, 3)
    
    dist.maha(x)
    #         1        2        3        4        5
    #2 2.362109                                    
    #3 1.725084 1.495655                           
    #4 2.959946 2.715641 2.690788                  
    #5 3.044610 1.218184 1.531026 2.717390         
    #6 2.740958 1.694767 2.877993 2.978265 2.794879
    

    结果与你的mhbd_calc2一致。

    【讨论】:

    • 那么,如果我理解正确的话,你 dist.maha 的精确度稍微低一些,但速度要快得多?精度为7位,和我的测试一样
    • 我可能错了,但 choleski 方法无法验证矩阵是否接近奇异。如果是这样,它可能会给出我们不想要的高值,不是吗?而solve() 进行此验证并返回错误以防止它。
    • 我认为这超出了我的知识范围,但我肯定会四处询问。另外,如果您不介意,您能否详细说明您的方法是如何工作的?这个功能肯定会为我节省很多时间,非常感谢:)
    • 嗨,我尝试通过使用 Optmatch 包中的 match_on 函数(默认是计算成对的马拉诺比斯距离)来交叉验证这个解决方案,我得到了不同的结果。你知道为什么会这样吗?
    猜你喜欢
    • 2015-06-18
    • 1970-01-01
    • 1970-01-01
    • 2016-04-21
    • 2013-09-10
    • 1970-01-01
    • 2023-03-12
    • 1970-01-01
    相关资源
    最近更新 更多