É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