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