Éste ejercicio es una reproducción del tutorial de Enhance Data Science para PCA con R. Los datos están tomados de la competencia en Kaggle de Digit Recognizer.
Leyendo los datos.
require(data.table)
mnist_data <- fread("C:/Users/Telemedicina/Documents/GENER/TESIS/data/MINST.csv")
Construyendo un grid de 3x3
par(mfrow = c(3,3))
for (i in 1:9)
{
##Changing i-th row to matrix
mat <- matrix(as.numeric(mnist_data[i,785:2, with = F]), nrow = 28, ncol = 28, byrow = F)
##Inverting row order
mat <- mat[nrow(mat):1,]
##plot
image(mat, main = paste0("Éste es un ", mnist_data[i,1,with=F]), col = grey(seq(1, 0, length = 256)))
}

Principal Components Analysis (PCA) para reducción de dimensionalidad de los datos.
Para el humano es relativamente sencillo reconocer patrones visuales, no obstante, desde el punto de vista de análisis de datos y procesamiento de los mismos, la cantidad de los mismos necesaria para representar la imagen anterior es significativa. Por ésta razón se vuelve relevante trabajar con la menor cantidad de datos posibles que nos sigan permietiendo obtener los resultados esperados, en el ejemplo anterior ésto sería el reconocer el número que está en la imagen. Para lograr ésto se utilizan técnicas de reducción de dimensionalidad de los datos, en éste caso el algoritmo PCA.
Preproceso de los datos.
Removing constant pixels.
El siguiente script muestra los pixeles que no tienen variabilidad en sus datos, como es de esperarse son los pixeles de los bordes. Por lo tanto hay que retirar los pixeles con valores constantes.
which(mnist_data[,sapply(.SD, FUN = function(x){min(x)==max(x)}), .SDcols= 2:ncol(mnist_data)])
pixel0 pixel1 pixel2 pixel3 pixel4 pixel5 pixel6 pixel7 pixel8 pixel9 pixel10
1 2 3 4 5 6 7 8 9 10 11
pixel11 pixel16 pixel17 pixel18 pixel19 pixel20 pixel21 pixel22 pixel23 pixel24 pixel25
12 17 18 19 20 21 22 23 24 25 26
pixel26 pixel27 pixel28 pixel29 pixel30 pixel31 pixel52 pixel53 pixel54 pixel55 pixel56
27 28 29 30 31 32 53 54 55 56 57
pixel57 pixel82 pixel83 pixel84 pixel85 pixel111 pixel112 pixel139 pixel140 pixel141 pixel168
58 83 84 85 86 112 113 140 141 142 169
pixel196 pixel392 pixel420 pixel421 pixel448 pixel476 pixel532 pixel560 pixel644 pixel645 pixel671
197 393 421 422 449 477 533 561 645 646 672
pixel672 pixel673 pixel699 pixel700 pixel701 pixel727 pixel728 pixel729 pixel730 pixel731 pixel754
673 674 700 701 702 728 729 730 731 732 755
pixel755 pixel756 pixel757 pixel758 pixel759 pixel760 pixel780 pixel781 pixel782 pixel783
756 757 758 759 760 761 781 782 783 784
El siguiente paso será el remover las celdas que están localizadas más cerca de un rango de 2 celdas de los bordes. Ésto removerá todos los pixeles constantes sin afectar significativamente a los dígitos.
mnist_data <- mnist_data[,-(which((1:784)%%28<= 2|(1:784)%%28>=26|1:784%/%28<=2|1:784%/%28>=26)+1),with=F]
Implementando el Algoritmo
PCA1 <- prcomp(mnist_data[,(2:ncol(mnist_data)), with = F], center = T, scale. = F)
Once the calculation is done it is important to visualize the variances explained by the n amount of principal components.
plot(cumsum(PCA1$sdev)/sum(PCA1$sdev)*100,
main = "Proporción Cumulativa de la Varianza Explicada",
xlab = "Cantidad de Componentes Principales",
ylab = "% De Varianza Explicada")

Visualización del PCA
En primer lugar se debe generar una proyección de los datos en el espacio (dimensiones) creados por el PCA.
projected <- scale(mnist_data[,(2:ncol(mnist_data)), with=F], PCA1$center, PCA1$scale) %*%
PCA1$rotation
Regresando al espacio de pixeles con los datos generados por el PCA, utilizando los primeros 60 componentes principales:
##Keeping only three dimensions
n_dim <- 70
##Projecting the data back using only the 3 principal components
coord_x <- data.table(mnist_data$label, projected[,1:n_dim] %*%
t(PCA1$rotation)[1:n_dim,])
par(mfrow = c(7,7), mar = c(0.1,0.1,0.1,0.1))
##Plotting 36 observations
for (i in 1:49)
{
mat <- matrix(as.numeric(coord_x[i, 530:2,with = F]),
nrow = 23, ncol = 23, byrow = F)
mat <- mat[nrow(mat):1,]
image(mat, main = paste0("Éste es un ", coord_x[i,1,with = F]), col = grey(seq(1, 0, length = 256)))
}

LS0tDQp0aXRsZTogIlBDQSBlbiBkYXRvcyBNYXRyaWNpYWxlcyINCm91dHB1dDogDQogIGh0bWxfbm90ZWJvb2s6IA0KICAgIHRoZW1lOiByZWFkYWJsZQ0KICAgIHRvYzogeWVzDQphdXRob3I6IEdlbmVyIEF2aWzDqXMgUg0KLS0tDQoNCsOJc3RlIGVqZXJjaWNpbyBlcyB1bmEgcmVwcm9kdWNjacOzbiBkZWwgdHV0b3JpYWwgZGUgW0VuaGFuY2UgRGF0YSBTY2llbmNlIHBhcmEgUENBIGNvbiBSXShodHRwOi8vZW5oYW5jZWRhdGFzY2llbmNlLmNvbS8yMDE3LzA1LzA3L3ItYmFzaWNzLXBjYS1yLykuDQpMb3MgZGF0b3MgZXN0w6FuIHRvbWFkb3MgZGUgbGEgY29tcGV0ZW5jaWEgZW4gW0thZ2dsZV0oaHR0cHM6Ly93d3cua2FnZ2xlLmNvbS9jL2RpZ2l0LXJlY29nbml6ZXIvZGF0YSkgZGUgKkRpZ2l0IFJlY29nbml6ZXIqLg0KDQojI0xleWVuZG8gbG9zIGRhdG9zLg0KDQpgYGB7ciwgbWVzc2FnZT1GQUxTRSwgd2FybmluZz1GQUxTRX0NCnJlcXVpcmUoZGF0YS50YWJsZSkNCm1uaXN0X2RhdGEgPC0gZnJlYWQoIkM6L1VzZXJzL1RlbGVtZWRpY2luYS9Eb2N1bWVudHMvR0VORVIvVEVTSVMvZGF0YS9NSU5TVC5jc3YiKQ0KYGBgDQoNCiMjQ29uc3RydXllbmRvIHVuICpncmlkKiBkZSAzeDMNCg0KYGBge3J9DQpwYXIobWZyb3cgPSBjKDMsMykpDQpmb3IgKGkgaW4gMTo5KQ0Kew0KICAjI0NoYW5naW5nIGktdGggcm93IHRvIG1hdHJpeA0KICBtYXQgPC0gbWF0cml4KGFzLm51bWVyaWMobW5pc3RfZGF0YVtpLDc4NToyLCB3aXRoID0gRl0pLCBucm93ID0gMjgsIG5jb2wgPSAyOCwgYnlyb3cgPSBGKQ0KICANCiAgIyNJbnZlcnRpbmcgcm93IG9yZGVyDQogIG1hdCA8LSBtYXRbbnJvdyhtYXQpOjEsXQ0KICANCiAgIyNwbG90DQogIGltYWdlKG1hdCwgbWFpbiA9IHBhc3RlMCgiw4lzdGUgZXMgdW4gIiwgbW5pc3RfZGF0YVtpLDEsd2l0aD1GXSksIGNvbCA9IGdyZXkoc2VxKDEsIDAsIGxlbmd0aCA9IDI1NikpKQ0KfQ0KYGBgDQojIypQcmluY2lwYWwgQ29tcG9uZW50cyBBbmFseXNpcyAoUENBKSogcGFyYSByZWR1Y2Npw7NuIGRlIGRpbWVuc2lvbmFsaWRhZCBkZSBsb3MgZGF0b3MuDQoNClBhcmEgZWwgaHVtYW5vIGVzIHJlbGF0aXZhbWVudGUgc2VuY2lsbG8gcmVjb25vY2VyIHBhdHJvbmVzIHZpc3VhbGVzLCBubyBvYnN0YW50ZSwgZGVzZGUgZWwgcHVudG8gZGUgdmlzdGEgZGUgYW7DoWxpc2lzIGRlIGRhdG9zIHkgcHJvY2VzYW1pZW50byBkZSBsb3MgbWlzbW9zLCBsYSBjYW50aWRhZCBkZSBsb3MgbWlzbW9zIG5lY2VzYXJpYSBwYXJhIHJlcHJlc2VudGFyIGxhIGltYWdlbiBhbnRlcmlvciBlcyBzaWduaWZpY2F0aXZhLiBQb3Igw6lzdGEgcmF6w7NuIHNlIHZ1ZWx2ZSByZWxldmFudGUgdHJhYmFqYXIgY29uIGxhIG1lbm9yIGNhbnRpZGFkIGRlIGRhdG9zIHBvc2libGVzIHF1ZSBub3Mgc2lnYW4gcGVybWlldGllbmRvIG9idGVuZXIgbG9zIHJlc3VsdGFkb3MgZXNwZXJhZG9zLCBlbiBlbCBlamVtcGxvIGFudGVyaW9yIMOpc3RvIHNlcsOtYSBlbCByZWNvbm9jZXIgZWwgbsO6bWVybyBxdWUgZXN0w6EgZW4gbGEgaW1hZ2VuLiANClBhcmEgbG9ncmFyIMOpc3RvIHNlIHV0aWxpemFuIHTDqWNuaWNhcyBkZSByZWR1Y2Npw7NuIGRlIGRpbWVuc2lvbmFsaWRhZCBkZSBsb3MgZGF0b3MsIGVuIMOpc3RlIGNhc28gZWwgYWxnb3JpdG1vIFBDQS4NCg0KIyMjUHJlcHJvY2VzbyBkZSBsb3MgZGF0b3MuDQoNCiMjIyNSZW1vdmluZyBjb25zdGFudCBwaXhlbHMuDQoNCkVsIHNpZ3VpZW50ZSAqc2NyaXB0KiBtdWVzdHJhIGxvcyBwaXhlbGVzIHF1ZSBubyB0aWVuZW4gdmFyaWFiaWxpZGFkIGVuIHN1cyBkYXRvcywgY29tbyBlcyBkZSBlc3BlcmFyc2Ugc29uIGxvcyBwaXhlbGVzIGRlIGxvcyBib3JkZXMuIFBvciBsbyB0YW50byBoYXkgcXVlIHJldGlyYXIgbG9zIHBpeGVsZXMgY29uIHZhbG9yZXMgY29uc3RhbnRlcy4NCmBgYHtyLCBlY2hvPVRSVUUsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9DQp3aGljaChtbmlzdF9kYXRhWyxzYXBwbHkoLlNELCBGVU4gPSBmdW5jdGlvbih4KXttaW4oeCk9PW1heCh4KX0pLCAuU0Rjb2xzPSAyOm5jb2wobW5pc3RfZGF0YSldKQ0KYGBgDQpFbCBzaWd1aWVudGUgcGFzbyBzZXLDoSBlbCByZW1vdmVyIGxhcyBjZWxkYXMgcXVlIGVzdMOhbiBsb2NhbGl6YWRhcyBtw6FzIGNlcmNhIGRlIHVuIHJhbmdvIGRlIDIgY2VsZGFzIGRlIGxvcyBib3JkZXMuIMOJc3RvIHJlbW92ZXLDoSB0b2RvcyBsb3MgcGl4ZWxlcyBjb25zdGFudGVzIHNpbiBhZmVjdGFyIHNpZ25pZmljYXRpdmFtZW50ZSBhIGxvcyBkw61naXRvcy4NCmBgYHtyLCBlY2hvPVRSVUUsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9DQptbmlzdF9kYXRhIDwtIG1uaXN0X2RhdGFbLC0od2hpY2goKDE6Nzg0KSUlMjg8PSAyfCgxOjc4NCklJTI4Pj0yNnwxOjc4NCUvJTI4PD0yfDE6Nzg0JS8lMjg+PTI2KSsxKSx3aXRoPUZdDQpgYGANCg0KIyMjSW1wbGVtZW50YW5kbyBlbCBBbGdvcml0bW8NCg0KYGBge3IsIG1lc3NhZ2U9RkFMU0UsIHdhcm5pbmc9RkFMU0V9DQpQQ0ExIDwtIHByY29tcChtbmlzdF9kYXRhWywoMjpuY29sKG1uaXN0X2RhdGEpKSwgd2l0aCA9IEZdLCBjZW50ZXIgPSBULCBzY2FsZS4gPSBGKQ0KYGBgDQpPbmNlIHRoZSBjYWxjdWxhdGlvbiBpcyBkb25lIGl0IGlzIGltcG9ydGFudCB0byB2aXN1YWxpemUgdGhlIHZhcmlhbmNlcyBleHBsYWluZWQgYnkgdGhlICpuKiBhbW91bnQgb2YgcHJpbmNpcGFsIGNvbXBvbmVudHMuDQpgYGB7cn0NCnBsb3QoY3Vtc3VtKFBDQTEkc2Rldikvc3VtKFBDQTEkc2RldikqMTAwLA0KICAgICBtYWluID0gIlByb3BvcmNpw7NuIEN1bXVsYXRpdmEgZGUgbGEgVmFyaWFuemEgRXhwbGljYWRhIiwNCiAgICAgeGxhYiA9ICJDYW50aWRhZCBkZSBDb21wb25lbnRlcyBQcmluY2lwYWxlcyIsDQogICAgIHlsYWIgPSAiJSBEZSBWYXJpYW56YSBFeHBsaWNhZGEiKQ0KYGBgDQoNCiMjI1Zpc3VhbGl6YWNpw7NuIGRlbCBQQ0ENCkVuIHByaW1lciBsdWdhciBzZSBkZWJlIGdlbmVyYXIgdW5hIHByb3llY2Npw7NuIGRlIGxvcyBkYXRvcyBlbiBlbCBlc3BhY2lvIChkaW1lbnNpb25lcykgY3JlYWRvcyBwb3IgZWwgUENBLg0KYGBge3J9DQpwcm9qZWN0ZWQgPC0gc2NhbGUobW5pc3RfZGF0YVssKDI6bmNvbChtbmlzdF9kYXRhKSksIHdpdGg9Rl0sIFBDQTEkY2VudGVyLCBQQ0ExJHNjYWxlKSAlKiUNCiAgUENBMSRyb3RhdGlvbg0KYGBgDQpSZWdyZXNhbmRvIGFsIGVzcGFjaW8gZGUgcGl4ZWxlcyBjb24gbG9zIGRhdG9zIGdlbmVyYWRvcyBwb3IgZWwgUENBLCB1dGlsaXphbmRvIGxvcyBwcmltZXJvcyA2MCBjb21wb25lbnRlcyBwcmluY2lwYWxlczoNCmBgYHtyfQ0KIyNLZWVwaW5nIG9ubHkgdGhyZWUgZGltZW5zaW9ucw0Kbl9kaW0gPC0gNzANCg0KIyNQcm9qZWN0aW5nIHRoZSBkYXRhIGJhY2sgdXNpbmcgb25seSB0aGUgMyBwcmluY2lwYWwgY29tcG9uZW50cw0KY29vcmRfeCA8LSBkYXRhLnRhYmxlKG1uaXN0X2RhdGEkbGFiZWwsIHByb2plY3RlZFssMTpuX2RpbV0gJSolDQogICAgICAgICAgICAgICAgICAgICAgICB0KFBDQTEkcm90YXRpb24pWzE6bl9kaW0sXSkNCnBhcihtZnJvdyA9IGMoNyw3KSwgbWFyID0gYygwLjEsMC4xLDAuMSwwLjEpKQ0KDQojI1Bsb3R0aW5nIDM2IG9ic2VydmF0aW9ucw0KZm9yIChpIGluIDE6NDkpDQp7DQogIG1hdCA8LSBtYXRyaXgoYXMubnVtZXJpYyhjb29yZF94W2ksIDUzMDoyLHdpdGggPSBGXSksDQogICAgICAgICAgICAgICAgbnJvdyA9IDIzLCBuY29sID0gMjMsIGJ5cm93ID0gRikNCiAgbWF0IDwtIG1hdFtucm93KG1hdCk6MSxdDQogIGltYWdlKG1hdCwgbWFpbiA9IHBhc3RlMCgiw4lzdGUgZXMgdW4gIiwgY29vcmRfeFtpLDEsd2l0aCA9IEZdKSwgY29sID0gZ3JleShzZXEoMSwgMCwgbGVuZ3RoID0gMjU2KSkpDQp9DQpgYGANCg0K