2016-12-07 41 views
1

我想要計算每個觀測值之間的數據集dat之間的Mahalanobis距離,其中每一行是一個觀測值,每一列都是一個變量。這樣的距離定義爲:每對觀測值的馬氏距離

formula

我寫的,做它的功能,但我覺得它是緩慢的。有沒有更好的方法來計算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秒)

回答

2

你需要一些數學知識。

  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同意。

+0

所以,如果我理解正確,你dist.maha稍微不夠精確,但更快?精度爲7位,與我的測試相同 – Oligg

+0

我可能是錯的,但是choleski方法不能驗證矩陣是否幾乎是單數。如果是這樣,它可以給我們不想要的高價值,不是嗎?而solve()會執行此驗證並返回一個錯誤以防止它。 – Oligg

+0

我認爲這超出了我的知識範圍,但我肯定會問。另外,如果你不介意,你可否詳細說明你的方法是如何工作的?這個功能肯定會節省我很多時間,非常感謝:) – Oligg