Implémentation étape par étape de PCA dans R à l'aide du didacticiel de Lindsay Smith

Je travaille dans R grâce à un excellent tutoriel PCA par Lindsay I Smith et je suis coincé dans la dernière étape. Le script R ci-dessous nous amène à l'étape (à la p.19) où les données originales sont reconstruites à partir de la composante principale (singulière dans ce cas), ce qui devrait produire un tracé en ligne droite le long de l'axe PCA1 (étant donné que les données n'a que 2 dimensions, dont la seconde est intentionnellement supprimée).

d = data.frame(x=c(2.5,0.5,2.2,1.9,3.1,2.3,2.0,1.0,1.5,1.1),
               y=c(2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9))

# mean-adjusted values 
d$x_adj = d$x - mean(d$x)
d$y_adj = d$y - mean(d$y)

# calculate covariance matrix and eigenvectors/values
(cm = cov(d[,1:2]))

#### outputs #############
#          x         y
# x 0.6165556 0.6154444
# y 0.6154444 0.7165556
##########################

(e = eigen(cm))

##### outputs ##############
# $values
# [1] 1.2840277 0.0490834
#
# $vectors
#          [,1]       [,2]
# [1,] 0.6778734 -0.7351787
# [2,] 0.7351787  0.6778734
###########################


# principal component vector slopes
s1 = e$vectors[1,1] / e$vectors[2,1] # PC1
s2 = e$vectors[1,2] / e$vectors[2,2] # PC2

plot(d$x_adj, d$y_adj, asp=T, pch=16, xlab='x', ylab='y')
abline(a=0, b=s1, col='red')
abline(a=0, b=s2)

# PCA data = rowFeatureVector (transposed eigenvectors) * RowDataAdjust (mean adjusted, also transposed)
feat_vec = t(e$vectors)
row_data_adj = t(d[,3:4])
final_data = data.frame(t(feat_vec %*% row_data_adj)) # ?matmult for details
names(final_data) = c('x','y')

#### outputs ###############
# final_data
#              x           y
# 1   0.82797019 -0.17511531
# 2  -1.77758033  0.14285723
# 3   0.99219749  0.38437499
# 4   0.27421042  0.13041721
# 5   1.67580142 -0.20949846
# 6   0.91294910  0.17528244
# 7  -0.09910944 -0.34982470
# 8  -1.14457216  0.04641726
# 9  -0.43804614  0.01776463
# 10 -1.22382056 -0.16267529
############################

# final_data[[1]] = -final_data[[1]] # for some reason the x-axis data is negative the tutorial's result

plot(final_data, asp=T, xlab='PCA 1', ylab='PCA 2', pch=16)

C'est aussi loin que j'ai, et tout va bien jusqu'à présent. Mais je ne peux pas comprendre comment les données sont obtenues pour le tracé final - la variance attribuable à l'ACP 1 - que Smith trace comme:

Voici ce que j'ai essayé (qui ignore l'ajout des moyens d'origine):

trans_data = final_data
trans_data[,2] = 0
row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))
plot(row_orig_data, asp=T, pch=16)

.. et a obtenu une erreur:

.. parce que j'ai perdu une dimension de données d'une manière ou d'une autre dans la multiplication matricielle. Je serais très reconnaissant d'avoir une idée de ce qui ne va pas ici.

* Éditer *

Je me demande si c'est la bonne formule:

row_orig_data = t(t(feat_vec) %*% t(trans_data))
plot(row_orig_data, asp=T, pch=16, cex=.5)
abline(a=0, b=s1, col='red')

Mais je suis un peu confus si oui car (a) je comprends les rowVectorFeaturebesoins à réduire à la dimensionnalité souhaitée (le vecteur propre pour PCA1), et (b) il ne correspond pas à l'abline PCA1:

Toutes vues très appréciées.

Vous étiez très près de là et vous avez été pris par un problème subtil en travaillant avec des matrices dans R. J'ai travaillé avec vous final_dataet j'ai obtenu les bons résultats indépendamment. Ensuite, j'ai regardé de plus près votre code. Pour abréger une longue histoire, où vous avez écrit

row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))

tu aurais été d'accord si tu avais écrit

row_orig_data = t(t(feat_vec) %*% t(trans_data))

trans_data$2\times1$$2\times10$t(feat_vec[1,])$1\times2$row_orig_data = t(as.matrix(feat_vec[1,],ncol=1,nrow=2) %*% t(trans_data))non-conformable arguments

row_orig_data = t(as.matrix(feat_vec[1,],ncol=1,nrow=2) %*% t(trans_data)[1,])

$2\times1$$1\times10$final_data$20=2\times10$row_orig_data$12=2\times1 + 1\times10$

$(XY)^T=Y^TX^T$t(t(p) %*% t(q)) = q %*% t

$x/y$$y/x$

Écrire

d_in_new_basis = as.matrix(final_data)

puis pour récupérer vos données dans leur base d'origine, vous avez besoin

d_in_original_basis = d_in_new_basis %*% feat_vec

Vous pouvez mettre à zéro les parties de vos données qui sont projetées le long du deuxième composant à l'aide de

d_in_new_basis_approx = d_in_new_basis
d_in_new_basis_approx[,2] = 0

et vous pouvez ensuite transformer comme avant

d_in_original_basis_approx = d_in_new_basis_approx %*% feat_vec

Le fait de les tracer sur le même tracé, ainsi que la ligne des composants principaux en vert, vous montre comment fonctionne l'approximation.

plot(x=d_in_original_basis[,1]+mean(d$x),
     y=d_in_original_basis[,2]+mean(d$y),
     pch=16, xlab="x", ylab="y", xlim=c(0,3.5),ylim=c(0,3.5),
     main="black=original data\nred=original data restored using only a single eigenvector")
points(x=d_in_original_basis_approx[,1]+mean(d$x),
       y=d_in_original_basis_approx[,2]+mean(d$y),
       pch=16,col="red")
points(x=c(mean(d$x)-e$vectors[1,1]*10,mean(d$x)+e$vectors[1,1]*10), c(y=mean(d$y)-e$vectors[2,1]*10,mean(d$y)+e$vectors[2,1]*10), type="l",col="green")

Revenons à ce que vous aviez. Cette ligne était ok

final_data = data.frame(t(feat_vec %*% row_data_adj))

feat_vec %*% row_data_adj$Y=S^TX$$S$$X$$Y$$Y$$X$$Y$$X$

Ensuite, vous aviez

trans_data = final_data
trans_data[,2] = 0

C'est correct: vous mettez à zéro les parties de vos données qui sont projetées le long du deuxième composant. Là où ça va mal, c'est

row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))

$\hat Y$$Y$$\mathbf{e}_1$t(feat_vec[1,]) %*% t(trans_data)$\mathbf{e}_1 \hat Y$

$2\times1$$2\times10$$\hat Y$$Y$$\mathbf{y}_1$$\mathbf{e}_1$$\mathbf{y}_1$$i$$\mathbf{e}_1\mathbf{y}_1$$\mathbf{e}_1$$i$

Je pense que vous avez la bonne idée, mais que vous êtes tombé sur une caractéristique désagréable de R. Ici encore, le morceau de code pertinent comme vous l'avez dit:

trans_data = final_data
trans_data[,2] = 0
row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))
plot(row_orig_data, asp=T, pch=16)

final_dataContient essentiellement les coordonnées des points d'origine par rapport au système de coordonnées défini par les vecteurs propres de la matrice de covariance. Pour reconstruire les points d'origine, il faut donc multiplier chaque vecteur propre avec la coordonnée transformée associée, par exemple

(1) final_data[1,1]*t(feat_vec[1,] + final_data[1,2]*t(feat_vec[2,])

ce qui donnerait les coordonnées d'origine du premier point. Dans votre question , vous définissez correctement le second composant à zéro, trans_data[,2] = 0. Si vous (comme vous l'avez déjà modifié) calculez

(2) row_orig_data = t(t(feat_vec) %*% t(trans_data))

vous calculez la formule (1) pour tous les points simultanément. Votre première approche

row_orig_data = t(t(feat_vec[1,]) %*% t(trans_data))

calcule quelque chose de différent et ne fonctionne que parce que R supprime automatiquement l'attribut dimension pour feat_vec[1,], il n'est donc plus un vecteur ligne mais traité comme un vecteur colonne. La transposition suivante en fait à nouveau un vecteur ligne et c'est la raison pour laquelle au moins le calcul ne produit pas d'erreur, mais si vous passez par le calcul, vous verrez que c'est quelque chose de différent de (1). En général, c'est une bonne idée dans les multiplications matricielles de supprimer la suppression de l'attribut de dimension qui peut être obtenue par le dropparamètre, par exemple feat_vec[1,,drop=FALSE].

$\Delta y / \Delta x$

s1 = e$vectors[2,1] / e$vectors[1,1] # PC1
s2 = e$vectors[2,2] / e$vectors[1,2] # PC2

Après avoir exploré cet exercice, vous pouvez essayer les méthodes les plus simples dans R. Il existe deux fonctions populaires pour effectuer l'ACP: princompet prcomp. La princompfonction effectue la décomposition des valeurs propres comme dans votre exercice. La prcompfonction utilise la décomposition en valeurs singulières. Les deux méthodes donneront les mêmes résultats presque tout le temps: cette réponse explique les différences de R, tandis que cette réponse explique les mathématiques . (Merci à TooTone pour les commentaires désormais intégrés à ce message.)

Ici, nous utilisons les deux pour reproduire l'exercice dans R. Tout d'abord en utilisant princomp:

d = data.frame(x=c(2.5,0.5,2.2,1.9,3.1,2.3,2.0,1.0,1.5,1.1), 
               y=c(2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9))

# compute PCs
p = princomp(d,center=TRUE,retx=TRUE)

# use loadings and scores to reproduce with only first PC
loadings = t(p$loadings[,1]) 
scores = p$scores[,1] 

reproduce = scores %*% loadings  + colMeans(d)

# plots
plot(reproduce,pch=3,ylim=c(-1,4),xlim=c(-1,4))
abline(h=0,v=0,lty=3)
mtext("Original data restored using only a single eigenvector",side=3,cex=0.7)

biplot(p)

Deuxième utilisation prcomp:

d = data.frame(x=c(2.5,0.5,2.2,1.9,3.1,2.3,2.0,1.0,1.5,1.1), 
               y=c(2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9))

# compute PCs
p = prcomp(d,center=TRUE,retx=TRUE)

# use loadings and scores to reproduce with only first PC
loadings = t(p$rotation[,1])
scores = p$x[,1]

reproduce = scores %*% loadings  + colMeans(d)

# plots
plot(reproduce,pch=3,ylim=c(-1,4),xlim=c(-1,4))
abline(h=0,v=0,lty=3)
mtext("Original data restored using only a single eigenvector",side=3,cex=0.7)

biplot(p)

De toute évidence, les signes sont inversés, mais l'explication de la variation est équivalente.