如何使用cor计算相关矩阵?

atmip9wb  于 2023-05-04  发布在  其他
关注(0)|答案(1)|浏览(176)

我用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位。

6xfqseft

6xfqseft1#

函数的第一行应该是这样的,因为R是逐列而不是逐行存储矩阵的

Q = t(t(X) - colMeans(X))

Q = X - matrix(colMeans(X), nrow(X), ncol(X), byrow = TRUE)

Q = scale(X, TRUE, FALSE)

甚至这个不是相同的Q,但最终给出了相同的答案

Q = scale(X)

如果我们从R的底部使用cov2cor,则

corr = function(X) {
  Q = scale(X)
  cov2cor(crossprod(Q))
}

相关问题