我在 R 中实现了自己的相关函数。令人惊讶的是,在使用内置
cor
函数时,我得到的结果略有不同。当 n
观察数量足够大时,差异似乎消失了。
我的职能:
corr = function(X) {
Q = X - colMeans(X)
S_ = colSums(Q**2)
S = sqrt(S_ %*% t(S_))
covarr = t(Q) %*% Q
corrr_ = covarr / S
return(corrr_)
}
library(mvtnorm)
set.seed(247)
X = rmvnorm(10, sigma = matrix(c(1,0.8,0.8,1), ncol=2)) # change 10 to 100, 1000, or 10000
corr(X)
cor(X)
对于
n=10
,我得到 0.8490966 与 0.8465363,所以变化在小数点后第三位。对于 n=1000
,我得到 0.7960206 与 0.7960925,所以变化在小数点后第 5 位。
函数的第一行应该是这样的,因为 R 是逐列而不是逐行存储矩阵
Q = t(t(X) - colMeans(X))
或
Q = scale(X, TRUE, FALSE)
甚至这个不是同一个Q但最终给出相同答案
Q = scale(X)
如果我们从 R 的基础上使用 cov2cor 那么
corr = function(X) {
Q = scale(X)
cov2cor(crossprod(Q))
}