Estudio sobre olivetti faces
El dataset de oliveti faces, consiste en un dataset de imágenes de caras tomadas en el año 1992 por el AT&T Laboratories en Cambridge. Hay 10 imágenes de cada uno de los sujetos, en diferentes ángulos, posiciones o sin con gafas. En total, tenemos este conjunto de 400 imágenes de dimensiones de 4096 (64x64). El objetivo del modelo en este dataset en comprobar si pertenece a una de estas 40 personas.
Este dataset fue creado con el objetivo de estudiar la capacidad de clasificar una imagen, en este caso una cara, usando el estudio de sus características principales. Este estudio, más tarde llamados Principal Component Analysis (PCA), son muy débiles cuando tienen interferencias de luz o ángulos, pero funcionan de manera excelente cuando las imágenes son uniformes, como es el caso de este dataset.
Dado que de cada persona sólo tenemos 10 imágenes tendremos que ver como actúa un modelo convolucional simple (el idóneo para imágenes). Por otro lado probaremos ciertas técnicas clásicas como Regresiones logísticas, SCM y estas mismas aplicando técnicas de extracción del componente principal (PCA).
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from sklearn.datasets import fetch_olivetti_faces
from sklearn.model_selection import train_test_split
from numpy.random import RandomState
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import classification_report
from sklearn.metrics import confusion_matrix
from sklearn.tree import DecisionTreeClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.model_selection import KFold
from sklearn.model_selection import cross_val_score
from sklearn.svm import SVC
from keras.utils import to_categorical
from sklearn.metrics import accuracy_score, precision_score, recall_score
from sklearn.model_selection import GridSearchCV
import seaborn
from sklearn.datasets import fetch_olivetti_faces
Conocer el dataset¶
Un primer paso para todo estudio es conocer dataset que tenemos, para ello en este caso es muy conveniente revisar qué datos tenemos, sus shapes, cuales son sus target y ver las imágenes con las que vamos a tratar.
n_row, n_col = 3, 10
n_components = n_row * n_col
image_shape = (64, 64)
def plot_gallery(title, images, n_col=n_col, n_row=n_row, cmap=plt.cm.gray):
plt.figure(figsize=(2. * n_col, 2.26 * n_row))
plt.suptitle(title, size=16)
for i, comp in enumerate(images):
plt.subplot(n_row, n_col, i + 1)
vmax = max(comp.max(), -comp.min())
plt.imshow(comp.reshape(image_shape), cmap=cmap,
interpolation='nearest',
vmin=-vmax, vmax=vmax)
plt.xticks(())
plt.yticks(())
plt.subplots_adjust(0.01, 0.05, 0.99, 0.93, 0.04, 0.)
#Formatear y validar
data = fetch_olivetti_faces()
X = data.data
X = X.reshape(X.shape[0], 64, 64, 1)
targets = data.target
print(targets)
#Veamos algunas personas del dataset
plot_gallery("First centered Olivetti faces", X[:n_components])
X_train, X_test, y_train, y_test=train_test_split(X, targets, test_size=0.2, stratify=targets, random_state=0)
print("X_train shape:",X_train.shape)
print("y_train shape:{}".format(y_train.shape))
#Veamos ahora algunas personas del dataset a entrenar
plot_gallery("First centered Olivetti faces", X_train[:n_components])
Creación del modelo¶
Vamos a empezar con unos modelos típicos.
import warnings
warnings.filterwarnings('ignore')
#Función para ver los resultados de un modelo simple
#Se puede customizar para que haga una validación cruzada
def show_result(model, x_train, y_train, y_train_pred, y_test, y_test_pred, do_cross):
print('Resultados en el conjunto de entrenamiento')
print(' Precisión:', accuracy_score(y_train, y_train_pred))
print(' Exactitud:', precision_score(y_train, y_train_pred, average = 'macro'))
print(' Exhaustividad:', recall_score(y_train, y_train_pred, average = 'macro'))
print('')
print(' Resultados en el conjunto de test')
print(' Precisión:', accuracy_score(y_test, y_test_pred))
print(' Exactitud:', precision_score(y_test, y_test_pred, average = 'macro'))
print(' Exhaustividad:', recall_score(y_test, y_test_pred, average = 'macro'))
if do_cross:
kfold=KFold(n_splits=5, shuffle=True, random_state=0)
cv_scores=cross_val_score(model, x_train, y_train, cv=kfold)
print("Cross validation score:{}".format(cv_scores.mean()))
cm = confusion_matrix(y_test, y_test_pred, labels=range(40))
plt.figure(figsize=(40, 40))
seaborn.heatmap(cm, annot = True, fmt = ".1f", linewidths = .5, cmap = 'Blues_r');
plt.ylabel('Actual label');
plt.xlabel('Predicted label');
Empezaremos con una regresión logística, la cual funciona realmente bien.
X_train_joined = X_train.reshape(X_train.shape[0],X_train.shape[1]*X_train.shape[2])
X_test_joined = X_test.reshape(X_test.shape[0],X_test.shape[1]*X_test.shape[2])
print("X_train shape:",X_train_joined.shape)
print("X_test shape:{}".format(X_test_joined.shape))
print(80*'-')
lr_mod = Pipeline([('std', StandardScaler()), ('lr', LogisticRegression(solver='lbfgs'))])
lr_mod.fit(X_train_joined, y_train)
print('LogisticRegression')
y_pred_trains = lr_mod.predict(X_train_joined)
y_pred_test = lr_mod.predict(X_test_joined)
show_result(lr_mod, X_train_joined, y_train, y_pred_trains, y_test, y_pred_test, True)
dtc_mod = DecisionTreeClassifier()
dtc_mod.fit(X_train_joined, y_train)
print('DecisionTreeClassifier')
y_pred_trains = dtc_mod.predict(X_train_joined)
y_pred_test = dtc_mod.predict(X_test_joined)
show_result(dtc_mod, X_train_joined, y_train, y_pred_trains, y_test, y_pred_test, False)
svc_mod = SVC(gamma='auto')
svc_mod.fit(X_train_joined, y_train)
print('SVM')
y_pred_trains = svc_mod.predict(X_train_joined)
y_pred_test = svc_mod.predict(X_test_joined)
show_result(svc_mod, X_train_joined, y_train, y_pred_trains, y_test, y_pred_test, False)
Bien hemos podido ver que con unas técnicas muy simple hemos obtenido unos buenos resultado, el que mejor ha funcionado ha sido la logistic regresion seguido de una SVM
Probemos ahora alguna técnica de búsqueda de transformación de variables o estudio de estas, como puede ser el caso de PCA
Principle Component Analysis (PCA) es un método que permite que los datos se representen en un tamaño menor. De acuerdo con este método, los datos se transforman en nuevos componentes y el tamaño de los datos se reduce al seleccionar los componentes más importantes.
Podemos estudiar la varianza de los datos para ver que número de componentes podemos elegir, esto también lo podemos hacer de forma paralela con el GridSearch CV. Aunque este último caso simplemente las buscaría por fuerza bruta.
from sklearn.decomposition import PCA
pca=PCA()
pca.fit(X_train_joined)
print(X_train_joined.shape)
plt.figure(1, figsize=(12,8))
plt.plot(pca.explained_variance_, linewidth=2)
plt.xlabel('Components')
plt.ylabel('Explained Variaces')
plt.show()
Probemos ahora para ver que resultados tenemos con PCA. Si no especificamos el número de componente este cogerá el mejor por defecto, luego veremos usando el GridSearch CV que este valor es 90, pero en la gráfica también se puede entrever. Ya que este es el punto en donde la recta tiende a aplanarse.
X_train_pca=pca.transform(X_train_joined)
X_test_pca=pca.transform(X_test_joined)
print(80*'-')
lr_mod = Pipeline([('lr', LogisticRegression())])
lr_mod.fit(X_train_pca, y_train)
print('LogisticRegression')
y_pred_trains = lr_mod.predict(X_train_pca)
y_pred_test = lr_mod.predict(X_test_pca)
show_result(lr_mod, X_train_pca, y_train, y_pred_trains, y_test, y_pred_test, False)
print(80*'-')
svc_mod = SVC(gamma='auto')
svc_mod.fit(X_train_pca, y_train)
print('SVM')
y_pred_trains = svc_mod.predict(X_train_pca)
y_pred_test = svc_mod.predict(X_test_pca)
show_result(svc_mod, X_train_pca, y_train, y_pred_trains, y_test, y_pred_test, False)
Como podemos ver, en ambos casos mejoramos en torno a un punto a los resultado sin PCA. Busquemos ahora los mejores resultados de PCA en número de componentes.
pipe = Pipeline(steps=[('pca', PCA()), ('logistic', LogisticRegression())])
param_grid = {
'pca__n_components': [5, 20, 30, 40, 50, 64, 90],
'logistic__penalty': ["l1","l2"] ,
}
estimator = GridSearchCV(pipe, param_grid)
estimator.fit(X_train_joined, y_train)
print("tuned hpyerparameters :(best parameters) ",estimator.best_params_)
Como podemos ver, confirmamos que el número de componentes es 90. Probemos ahora estos nuevos resultados.
pca=PCA(n_components=90)
pca.fit(X_train_joined)
X_train_pca=pca.transform(X_train_joined)
X_test_pca=pca.transform(X_test_joined)
print(80*'-')
lr_mod = Pipeline([('lr', LogisticRegression(penalty='l2'))])
lr_mod.fit(X_train_pca, y_train)
print('LogisticRegression')
y_pred_trains = lr_mod.predict(X_train_pca)
y_pred_test = lr_mod.predict(X_test_pca)
show_result(lr_mod, X_train_pca, y_train, y_pred_trains, y_test, y_pred_test, False)
print(80*'-')
svc_mod = SVC(gamma='auto')
svc_mod.fit(X_train_pca, y_train)
print('SVM')
y_pred_trains = svc_mod.predict(X_train_pca)
y_pred_test = svc_mod.predict(X_test_pca)
show_result(svc_mod, X_train_pca, y_train, y_pred_trains, y_test, y_pred_test, False)
Realmente, con unos datos de test de 0.96 de precisión y 0.979 con un simple regresión logistica + PCA son excelentes, y en un entorno de producción no haría falta seguir probando más. Ya que supondría invertir tiempo y dinero en una mejora muy pequeña, pero para este estudio vamos a probar por curiosidad con convolucionales. Estas están pensadas específicamente para imágenes, pero no funcionan tan bien con datasets pequeños.
Convolucionales¶
Vamos a probar ahora con una red profunda tenemos pocos datos, pero una red convolucional se podría adaptar bien, ya que estan creadas para esto. Como se puede ver hemos añadido una capa de dropout para evitar el overfitting.
from keras import layers
from keras import Model
from keras.optimizers import RMSprop
from keras.models import Sequential
from keras.layers import Dense, Conv2D, Flatten
#transfor data into input shape
X_train, X_test, y_train, y_test=train_test_split(X, targets, test_size=0.2, stratify=targets, random_state=0)
y_train = to_categorical(y_train)
y_test = to_categorical(y_test)
img_input = layers.Input(shape=(64, 64, 1))
# First convolution extracts 16 filters that are 3x3
# Convolution is followed by max-pooling layer with a 2x2 window
x = layers.Conv2D(16, 3, activation='relu')(img_input)
x = layers.MaxPooling2D(2)(x)
# Second convolution extracts 32 filters that are 3x3
# Convolution is followed by max-pooling layer with a 2x2 window
x = layers.Conv2D(32, 3, activation='relu')(x)
x = layers.MaxPooling2D(2)(x)
# Flatten feature map to a 1-dim tensor
x = layers.Flatten()(x)
# Create a fully connected layer with ReLU activation and 512 hidden units
x = layers.Dense(256, activation='relu')(x)
# Add a dropout rate of 0.5
x = layers.Dropout(0.2)(x)
# Create output layer with a single node and sigmoid activation
output = layers.Dense(40, activation='softmax')(x)
# Configure and compile the model
model = Model(img_input, output)
model.summary()
Train del modelo¶
from keras import optimizers
rms = optimizers.RMSprop(lr=0.001, rho=0.9, epsilon=None, decay=0.0)
model.compile(loss='categorical_crossentropy',
optimizer=rms,
metrics=['mae','acc'])
history = model.fit(X_train, y_train, validation_data=(X_test, y_test), epochs=25)
Evaluación del modelo¶
from sklearn.metrics import precision_recall_curve, roc_curve, auc
# Retrieve a list of accuracy results on training and test data
# sets for each training epoch
acc = history.history['acc']
val_acc = history.history['val_acc']
# Retrieve a list of list results on training and test data
# sets for each training epoch
loss = history.history['loss']
val_loss = history.history['val_loss']
# Get number of epochs
epochs = range(len(acc))
# Plot training and validation accuracy per epoch
plt.plot(epochs, acc)
plt.plot(epochs, val_acc)
plt.title('Training and validation accuracy')
plt.figure()
# Plot training and validation loss per epoch
plt.plot(epochs, loss)
plt.plot(epochs, val_loss)
y_pred_keras = model.predict(X_test).ravel()
y_pred_keras = y_pred_keras.reshape(y_test.shape)
# Compute ROC curve and ROC area for each class
fpr = dict()
tpr = dict()
roc_auc = dict()
for i in range(40):
fpr[i], tpr[i], _ = roc_curve(y_test[:, i], y_pred_keras[:, i])
roc_auc[i] = auc(fpr[i], tpr[i])
# Compute micro-average ROC curve and ROC area
fpr["micro"], tpr["micro"], _ = roc_curve(y_test.ravel(), y_pred_keras.ravel())
roc_auc["micro"] = auc(fpr["micro"], tpr["micro"])
from scipy import interp
plt.figure()
lw = 2
plt.plot(fpr[4], tpr[4], color='darkorange',
lw=lw, label='ROC curve (area = %0.2f)' % roc_auc[2])
plt.plot([0, 1], [0, 1], color='navy', lw=lw, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Receiver operating characteristic example')
plt.legend(loc="lower right")
plt.show()
n_classes = 40
# First aggregate all false positive rates
all_fpr = np.unique(np.concatenate([fpr[i] for i in range(n_classes)]))
# Then interpolate all ROC curves at this points
mean_tpr = np.zeros_like(all_fpr)
for i in range(n_classes):
mean_tpr += interp(all_fpr, fpr[i], tpr[i])
# Finally average it and compute AUC
mean_tpr /= n_classes
fpr["macro"] = all_fpr
tpr["macro"] = mean_tpr
roc_auc["macro"] = auc(fpr["macro"], tpr["macro"])
plt.figure()
plt.plot(fpr["micro"], tpr["micro"],
label='micro-average ROC curve (area = {0:0.2f})'
''.format(roc_auc["micro"]),
color='red', linewidth=4)
plt.plot(fpr["macro"], tpr["macro"],
label='macro-average ROC curve (area = {0:0.2f})'
''.format(roc_auc["macro"]),
color='blue', linewidth=4)
plt.plot([0, 1], [0, 1], 'k--', lw=lw)
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Some extension of Receiver operating characteristic to multi-class')
plt.legend(loc="lower right")
plt.show()
plt.figure(2)
plt.xlim(0, 0.2)
plt.ylim(0.8, 1)
plt.plot([0, 1], [0, 1], 'k--')
plt.plot(fpr["micro"], tpr["micro"],label='micro-average ROC curve (area = {0:0.2f})'.format(roc_auc["micro"]),color='red', linewidth=4)
plt.plot(fpr["macro"], tpr["macro"], label='macro-average ROC curve (area = {0:0.2f})'.format(roc_auc["macro"]), color='blue', linewidth=4)
plt.xlabel('False positive rate')
plt.ylabel('True positive rate')
plt.title('ROC curve (zoomed in at top left)')
plt.legend(loc='best')
plt.show()
score_train, mae, acc_train = model.evaluate(X_train, y_train)
score, mae, acc = model.evaluate(X_test, y_test)
print('Resultados en el conjunto de entrenamiento')
print(' Precisión:', acc_train)
print(' Score:', score_train)
print('')
print(' Resultados en el conjunto de test')
print(' Precisión:', acc)
print(' Score:', score)
Conclusiones¶
Las redes convolucionales tienen un futuro prometedor, estoy seguro que con data augmentation y mas epoch se mejoraría ese 0.95 de acc en test. Pero tal vez todo este esfuerzo no es necesario, con una simple linear regresion + PCA hemos obtenido el mejor resultado:
- Resultados en el conjunto de test
- Precisión: 0.9625
- Exactitud: 0.9791666666666667
- Exhaustividad: 0.9625
Muchas veces no es necesaria las técnicas más punteras, con un simples técnicas que llevan desde los años 90 podemos obtener excelentes resultados.
Comentar