2016-12-07 14 views
1

データセットdatの各観測間のマハラノビス距離を計算しようとしています。ここで、各行は観測値であり、各列は変数です。このような距離は次のように定義されます観測の各ペアのマハラノビス距離

formula

私はそれをしない機能を書いたが、それが遅いような気がします。これをRで計算する方法はありますか?

は、機能をテストするために、いくつかのデータを生成するには:

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

これは私がこれまでに書かれている機能です。これらは両方とも機能し、データ(800個の変数と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)) 
} 

NB:私はouter()を使用してみましたが、それはさらに遅い(60秒)だった

答えて

2

あなたは、いくつかの数学の知識が必要です。

  1. 経験的共分散のコレスキー分解を行い、観測値を標準化します。
  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

私は間違っているかもしれませんが、コレスキーメソッドは行列がほぼ特異であるかどうかを検証しません。それが事実なら、私たちが望んでいない高い価値を与えるかもしれない、いいえ? solve()はこの検証を行い、それを防ぐためにエラーを返します。 – Oligg

+0

それは私の知識の限界を超えていると思うが、確かに尋ねるだろう。また、あなたが気にしないであなたの方法がどのように機能するかについて少し詳しく説明できますか?この機能は、確かに大きな時間を節約します、ありがとうalot :) – Oligg

関連する問題