¿Paso por aplicación de paso de PCA en R: Cómo reconstruir los datos originales del primera componente principal?

Question

Estoy trabajando en R a través de una excelente PCA tutorial de Lindsay me Smith y estoy estancado en la última etapa. La R siguiente script nos lleva hasta el escenario (en la p.19), en donde los datos originales han sido reconstruidos a partir de la (singular, en este caso) de la Componente Principal, que debe producir una línea recta de la trama a lo largo de la PCA1 eje (dado que los datos sólo tiene 2 dimensiones, el segundo de los cuales está siendo intencionalmente ha caído).

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)

Esto es lo que he conseguido, y todo OK hasta el momento. Pero no puedo averiguar cómo se obtienen los datos para la final de la trama - la varianza atribuible a la PCA 1 - Smith parcelas como:

Esto es lo que he intentado (que ignora la adición al medio original):

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)

.. y tengo un erronous:

.. porque he perdido a una dimensión de datos de alguna manera en la multiplicación de la matriz. Yo estaría muy agradecido por una idea de lo que va mal aquí.

---

* Edit *

Me pregunto si esta es la fórmula correcta:

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

Pero estoy un poco confundido si es así debido a que (a) entiendo que el rowVectorFeature debe ser reducida a la deseada dimensionalidad (el vector propio para PCA1), y (b) no se alinean con la PCA1 abline:

Cualquier punto de vista muy apreciado.

Answer 1

Que eran muy, muy cerca de allí y quedó atrapado por un problema sutil en el trabajo con matrices en R. he trabajado a través de su final_data y obtuvo los resultados correctos de forma independiente. Entonces yo tenía una mirada más cercana a su código. Para cortar una larga historia corta, en la que escribió

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

usted hubiera estado bien si hubiera escrito

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

en lugar de ello (porque iba a ceros la parte de trans_data que se proyecta sobre el segundo vector propio). Como era usted tratando de multiplicar un $2\times1$ matriz por un $2\times10$ matrix pero R no dará un error. El problema es que t(feat_vec[1,]) es tratado como $1\times2$. Tratando row_orig_data = t(as.matrix(feat_vec[1,],ncol=1,nrow=2) %*% t(trans_data)) hubiera dado un non-conformable arguments de error. El siguiente, posiblemente más a lo largo de las líneas de lo que pretende, también han trabajado

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

como se multiplica un $2\times1$ matriz por un $1\times10$ matriz (tenga en cuenta que usted podría haber usado el original final_data de la matriz de aquí). No es necesario hacerlo de esta manera, pero es más agradable matemáticamente porque demuestra que usted está recibiendo $20=2\times10$ valores en row_orig_data de $12=2\times1 + 1\times10$ valores en el lado derecho.

He dejado mi respuesta original a continuación, como alguien podría encontrar que es útil, y que demuestra conseguir la necesaria parcelas. También muestra que el código puede ser un poco más simple por deshacerse de algunos innecesarios transpone: $(XY)^T=Y^TX^T$ t(t(p) %*% t(q)) = q %*% t.

Volver a la edición, he añadido el componente principal de la línea en verde para mi parcela a continuación. En su pregunta se obtuvo tenía la pendiente en$x/y$$y/x$.

---

Escribir

d_in_new_basis = as.matrix(final_data)

a continuación, para recuperar los datos en su base original que usted necesita

d_in_original_basis = d_in_new_basis %*% feat_vec

Usted puede poner en cero las partes de los datos que se proyecta a lo largo de la segunda componente mediante

d_in_new_basis_approx = d_in_new_basis
d_in_new_basis_approx[,2] = 0

y entonces usted puede transformar como antes

d_in_original_basis_approx = d_in_new_basis_approx %*% feat_vec

Trazado de estas en la misma parcela, junto con el componente principal de la línea en verde, se muestra cómo la aproximación trabajado.

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

Recapitulemos lo que tenía. Esta línea fue ok

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

El crucial poco aquí es feat_vec %*% row_data_adj que es equivalente a $Y=S^TX$ donde $S$ es la matriz de vectores propios y $X$ es la matriz de datos con los datos en filas y $Y$ es los datos en la nueva base. Lo que esto dice es que la primera fila de $Y$ es la suma de (filas de $X$ ponderado por el primer autovector). Y la segunda fila de $Y$ es la suma de (filas de $X$ ponderado por el segundo vector propio).

Entonces tenía

trans_data = final_data
trans_data[,2] = 0

Esto está bien: eres sólo la reducción a cero de las partes de los datos que se proyecta a lo largo de la segunda componente. Donde las cosas van mal es

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

Escrito $\hat Y$ para la matriz de datos $Y$ en la nueva base, con ceros en la segunda fila, y escribir $\mathbf{e}_1$ para el primer vector propio, el final del negocio de este código t(feat_vec[1,]) %*% t(trans_data) trata $\mathbf{e}_1 \hat Y$.

Como se explicó anteriormente (aquí es donde me di cuenta de la sutil R problema y escribió la primera parte de mi respuesta), matemáticamente usted está tratando de multiplicar un $2\times1$ vector por un $2\times10$ matriz. Esto no funciona matemáticamente. Lo que usted debe hacer es tomar la primera fila de $\hat Y$ = la primera fila de $Y$: llame a este $\mathbf{y}_1$. Después multiplica $\mathbf{e}_1$ $\mathbf{y}_1$ juntos. El $i$ésima columna del resultado $\mathbf{e}_1\mathbf{y}_1$ es el autovector $\mathbf{e}_1$ ponderado por el 1 de coordenadas sólo de la $i$th punto en la nueva base, que es lo que quieres.

Answer 2

Creo que tienes la idea correcta, pero se tropezó con una desagradable característica de R. Aquí de nuevo el código correspondiente pieza, como se ha declarado:

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)

Esencialmente final_data contiene las coordenadas de los puntos originales con respecto al sistema de coordenadas definido por los vectores propios de la matriz de covarianza. Para reconstruir los puntos originales por lo tanto, se ha de multiplicar cada autovector asociado transformadas de coordenadas, por ejemplo,

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

que llevaría a la original coordenadas del primer punto. En su pregunta de establecer el segundo componente correctamente a cero, trans_data[,2] = 0. Si entonces (como ya editado) calcular

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

calcular la fórmula (1) para todos los puntos simultáneamente. Su primer acercamiento

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

calcula algo diferente y sólo funciona porque R automáticamente cae el atributo de la dimensión de feat_vec[1,], por lo que no es un vector fila más, pero tratada como un vector columna. La posterior transposición convierte en un vector fila de nuevo y esa es la razón por la que, al menos, el cálculo no se producirá ningún error, pero si usted va a través de las matemáticas que se ve que es algo diferente (1). En general es una buena idea en la matriz de multiplicaciones para suprimir la caída de el atributo de la dimensión que pueden ser alcanzados por la drop parámetro, por ejemplo feat_vec[1,,drop=FALSE].

Tu editado solución parece la correcta, pero calcula la pendiente, si PCA1 erróneamente. La pendiente está dada por $\Delta y / \Delta x$, por lo tanto

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

Answer 3

Después de explorar este ejercicio, usted puede tratar de las formas más fáciles en R. Hay dos funciones populares para hacer PCA: princomp y prcomp. El princomp función hace que el autovalor de la descomposición como se ha hecho en su ejercicio. El prcomp función utiliza la descomposición de valor singular. Ambos métodos dan el mismo resultado casi todo el tiempo: esta respuesta explica las diferencias en R, mientras que esta respuesta explica las matemáticas. (Gracias a TooTone para comentarios, ahora integrado en este post.)

Aquí se utiliza tanto para reproducir el ejercicio en R. en Primer lugar usando 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)

Segundo, con 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)

Claramente los signos se da la vuelta, pero la explicación de la variación equivalente.