自行计算协方差矩阵（不使用 `cov`）

Question

我正在学习有关协方差矩阵的教程，可以在此处找到：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) %[=27 究竟是多少=]%D正在积累。

有人能给我解释一下吗？

Answer 1

But I do not get how exactly t(D) %% D is built up.

这是矩阵叉积，矩阵乘法的一种特殊形式。如果你不明白它在做什么，考虑下面的 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

自行计算协方差矩阵（不使用 `cov`）

Compute covariance matrix on our own (without using `cov`)

r

matrix

covariance