我正在遵循一个关于协方差矩阵的教程,可以在这里找到:http://stats.seandolinar.com/making-a-covariance-matrix-in-r/
包括以下步骤:
#create a dataframe
a <- c(1,2,3,4,5,6)
b <- c(2,3,5,6,1,9)
c <- c(3,5,5,5,10,8)
d <- c(10,20,30,40,50,55)
e <- c(7,8,9,4,6,10)
#create matrix from vectors
M <- cbind(a,b,c,d,e)
M_mean <- matrix(data=1, nrow=n) %*% cbind(mean(a),mean(b),mean(c),mean(d),mean(e))
k <- ncol(M) #number of variables
n <- nrow(M) #number of subjects
然后创建一个像这样的差矩阵:
D <- M - M_mean
这对我来说很简单。但下一步是创建协方差矩阵:
C <- (n-1)^-1 t(D) %*% D
我得到t(D) %% D除以(n-1)^1 = 6。但是我不明白t(D) %% D是如何形成的。
有人能给我解释一下吗? 但是我不知道t(D) %% D究竟是如何形成的。
这是矩阵外积,矩阵乘法的一种特殊形式。如果你不明白它在做什么,考虑下面的R循环来帮助你理解它:
DtD <- matrix(0, nrow = ncol(D), ncol = ncol(D))
for (j in 1:ncol(D))
for (i in 1:ncol(D))
DtD[i, j] <- sum(D[, i] * D[, j])
注意,实际上没有人会为这个写R循环;这只是为了帮助你理解算法。
原始回答
假设我们有一个矩阵X
,其中每一列给出一个特定随机变量的观测值,通常我们只使用R基函数cov(X)
来得到协方差矩阵。
现在你想自己写一个协方差函数;这也不难(我很久以前就做过这个练习)。它需要3个步骤:
- 列居中(即对所有变量取均值);
- 矩阵外积;
- 平均(对
nrow(X) - 1
而不是nrow(X)
进行偏置调整)。
这个短代码做到了:
crossprod(sweep(X, 2L, colMeans(X))) / (nrow(X) - 1L)
考虑一个小示例
set.seed(0)
## 3 variable, each with 10 observations
X <- matrix(rnorm(30), nrow = 10, ncol = 3)
## reference computation by `cov`
cov(X)
# [,1] [,2] [,3]
#[1,] 1.4528358 -0.20093966 -0.10432388
#[2,] -0.2009397 0.46086672 -0.05828058
#[3,] -0.1043239 -0.05828058 0.48606879
## own implementation
crossprod(sweep(X, 2L, colMeans(X))) / (nrow(X) - 1L)
# [,1] [,2] [,3]
#[1,] 1.4528358 -0.20093966 -0.10432388
#[2,] -0.2009397 0.46086672 -0.05828058
#[3,] -0.1043239 -0.05828058 0.48606879
如果你想得到相关矩阵怎么办?
有很多方法。如果我们想直接获取它,执行:
crossprod(scale(X)) / (nrow(X) - 1L)
# [,1] [,2] [,3]
#[1,] 1.0000000 -0.2455668 -0.1241443
#[2,] -0.2455668 1.0000000 -0.1231367
#[3,] -0.1241443 -0.1231367 1.0000000
如果我们想首先获得协方差,然后(对称地)通过根对角线重新缩放它以获得相关性,我们可以这样做:
## covariance first
V <- crossprod(sweep(X, 2L, colMeans(X))) / (nrow(X) - 1L)
## symmetric rescaling
V / tcrossprod(diag(V) ^ 0.5)
# [,1] [,2] [,3]
#[1,] 1.0000000 -0.2455668 -0.1241443
#[2,] -0.2455668 1.0000000 -0.1231367
#[3,] -0.1241443 -0.1231367 1.0000000
我们还可以使用服务R函数cov2cor
将协方差转换为相关性:
cov2cor(V)
# [,1] [,2] [,3]
#[1,] 1.0000000 -0.2455668 -0.1241443
#[2,] -0.2455668 1.0000000 -0.1231367
#[3,] -0.1241443 -0.1231367 1.0000000