matrice de covariance Compute sur notre propre (sans utiliser `cov`)

Question

Je suivais un tutoriel sur les matrices de covariance qui pourraient être trouvés ici: http://stats.seandolinar.com/making-a-covariance-matrix-in-r/

Il comprend les étapes suivantes:

#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 

Et puis créer une matrice de différences comme celle-ci:

D <- M - M_mean 

Tout cela est très joli pour moi. Mais la prochaine étape que cela pour créer une matrice de covariance:

C <- (n-1)^-1 t(D) %*% D 

Je reçois que la partie t (D)% % D est divisé par (n-1)^1 = 6. Mais je ne suis pas comment exactement t (D)%% D est construit.

Quelqu'un pourrait-il m'expliquer cela?

Answer 1

Mais je ne comprends pas exactement comment t (D) %% D est construit.

Il s'agit d'un produit à matrice croisée, une forme spéciale de multiplication matricielle. Si vous ne comprenez pas ce qu'il fait, considérer la boucle R suivante pour vous aider à absorber ceci:

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]) 

Note, personne ne va en réalité écrire boucle R pour cela; C'est juste pour vous aider à comprendre l'algorithme.

---

originale Réponse

Supposons que nous ayons une matrice X, où chaque colonne donne des observations d'une variable aléatoire spécifique, normalement, nous utilisons simplement la fonction de base R cov(X) pour obtenir la matrice de covariance.

Maintenant, vous voulez écrire une fonction de covariance vous-même; Ce n'est pas difficile non plus (je l'ai fait il y a longtemps comme exercice). Il faut 3 étapes:

• centrage de colonne (c'est-à-dire, dé-moyenne pour toutes les variables);
• produit croisé de matrice;
• moyenne (sur nrow(X) - 1 et non nrow(X) pour l'ajustement de la polarisation).

Ce code court le fait:

crossprod(sweep(X, 2L, colMeans(X)))/(nrow(X) - 1L) 

Tenir compte un petit exemple

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 

---

Que faire si vous voulez obtenir la matrice de corrélation?

Il existe plusieurs façons.Si nous voulons obtenir directement, faire:

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 

Si nous voulons d'abord obtenir covariance, puis (symétriquement) RESCALE par la racine diagonale pour obtenir une corrélation, nous pouvons faire:

## 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 

Nous peut également utiliser une fonction de service de R cov2cor pour convertir covariance à la corrélation:

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