Yo creo que lo que usted está consiguiendo en su pregunta se refiere al truncamiento de los datos utilizando un número menor de componentes principales (PC). Para este tipo de operaciones, creo que la función prcomp
es más ilustrativa en la que es más fácil visualizar la multiplicación de la matriz utilizada en la reconstrucción.
En primer lugar, una síntesis de conjunto de datos, Xt
, de llevar a cabo la PCA (normalmente, usted tendría que centrar las muestras con el fin de describir PC relativas a una matriz de covarianza:
#Generate data
m=50
n=100
frac.gaps <- 0.5 # the fraction of data with NaNs
N.S.ratio <- 0.25 # the Noise to Signal ratio for adding noise to data
x <- (seq(m)*2*pi)/m
t <- (seq(n)*2*pi)/n
#True field
Xt <-
outer(sin(x), sin(t)) +
outer(sin(2.1*x), sin(2.1*t)) +
outer(sin(3.1*x), sin(3.1*t)) +
outer(tanh(x), cos(t)) +
outer(tanh(2*x), cos(2.1*t)) +
outer(tanh(4*x), cos(0.1*t)) +
outer(tanh(2.4*x), cos(1.1*t)) +
tanh(outer(x, t, FUN="+")) +
tanh(outer(x, 2*t, FUN="+"))
Xt <- t(Xt)
#PCA
res <- prcomp(Xt, center = TRUE, scale = FALSE)
names(res)
En los resultados o prcomp
, se puede ver la PC (res$x
), los valores propios (res$sdev
), dando información sobre la magnitud de cada PC, y las cargas (res$rotation
).
res$sdev
length(res$sdev)
res$rotation
dim(res$rotation)
res$x
dim(res$x)
Elevando al cuadrado los valores propios, se obtiene la varianza explicada por cada PC:
plot(cumsum(res$sdev^2/sum(res$sdev^2))) #cumulative explained variance
Por último, puede crear una versión truncada de sus datos mediante el uso de sólo el líder (importante) PCs:
pc.use <- 3 # explains 93% of variance
trunc <- res$x[,1:pc.use] %*% t(res$rotation[,1:pc.use])
#and add the center (and re-scale) back to data
if(res$scale != FALSE){
trunc <- scale(trunc, center = FALSE , scale=1/res$scale)
}
if(res$center != FALSE){
trunc <- scale(trunc, center = -1 * res$center, scale=FALSE)
}
dim(trunc); dim(Xt)
Se puede ver que el resultado es un poco más suave matriz de datos, a pequeña escala las características de filtrado:
RAN <- range(cbind(Xt, trunc))
BREAKS <- seq(RAN[1], RAN[2],,100)
COLS <- rainbow(length(BREAKS)-1)
par(mfcol=c(1,2), mar=c(1,1,2,1))
image(Xt, main="Original matrix", xlab="", ylab="", xaxt="n", yaxt="n", breaks=BREAKS, col=COLS)
box()
image(trunc, main="Truncated matrix (3 PCs)", xlab="", ylab="", xaxt="n", yaxt="n", breaks=BREAKS, col=COLS)
box()
Y aquí es un enfoque básico que usted puede hacer fuera de la prcomp función:
#alternate approach
Xt.cen <- scale(Xt, center=TRUE, scale=FALSE)
C <- cov(Xt.cen, use="pair")
E <- svd(C)
A <- Xt.cen %*% E$u
#To remove units from principal components (A)
#function for the exponent of a matrix
"%^%" <- function(S, power)
with(eigen(S), vectors %*% (values^power * t(vectors)))
Asc <- A %*% (diag(E$d) %^% -0.5) # scaled principal components
#Relationship between eigenvalues from both approaches
plot(res$sdev^2, E$d) #PCA via a covariance matrix - the eigenvalues now hold variance, not stdev
abline(0,1) # same results
Ahora, decidir qué equipos a retener es otra cuestión - a la que yo estaba interesado en un tiempo. Espero que ayude.