1
我想要計算每個觀測值之間的數據集dat
之間的Mahalanobis距離,其中每一行是一個觀測值,每一列都是一個變量。這樣的距離定義爲:每對觀測值的馬氏距離
我寫的,做它的功能,但我覺得它是緩慢的。有沒有更好的方法來計算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秒)
所以,如果我理解正確,你dist.maha稍微不夠精確,但更快?精度爲7位,與我的測試相同 – Oligg
我可能是錯的,但是choleski方法不能驗證矩陣是否幾乎是單數。如果是這樣,它可以給我們不想要的高價值,不是嗎?而solve()會執行此驗證並返回一個錯誤以防止它。 – Oligg
我認爲這超出了我的知識範圍,但我肯定會問。另外,如果你不介意,你可否詳細說明你的方法是如何工作的?這個功能肯定會節省我很多時間,非常感謝:) – Oligg