introduction

Beaucoup de gens liront First Pattern Recognition comme introduction à l'apprentissage automatique. J'étais l'un d'entre eux, et je pensais l'avoir lu une fois et l'avoir compris complètement, mais quand je l'ai relu, je ne comprends toujours rien [^ 1], alors j'ai décidé de le relire en le mettant en œuvre. Cet article résume le code qui reproduit l'exemple d'exécution décrit dans le livre en Python (l'environnement d'exploitation du code décrit est numpy 1.16.2, pandas 0.24.2). Cette fois, c'est le chapitre 4, «Modèle de probabilité et fonction de discrimination».

Exemple d'exécution 4.1 Normalisation

La normalisation est la soustraction de la valeur moyenne et la division par l'écart type. Ci-dessous le code.

--Chargement des bibliothèques requises

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_style('whitegrid')

--Lecture des données d'iris

from sklearn.datasets import load_iris

iris = load_iris()
data = pd.DataFrame(iris['data'], columns=iris.feature_names)
target = pd.Series(iris['target']).map({0: 'setosa', 1:'versicolor', 2:'virginica'})

Longueur et largeur moyennes des pétales d'iris et matrice de covariance

print('mu = ', np.mean(data[['petal length (cm)', 'petal width (cm)']]))
print('Sigma = ', np.cov(data[['petal length (cm)', 'petal width (cm)']], rowvar=0))

#production=>
# mu = 
#  petal length (cm)    3.758000
# petal width (cm)     1.199333
# dtype: float64
# Sigma = 
#  [[3.11627785 1.2956094 ]
#  [1.2956094  0.58100626]]

Diagramme de dispersion avant normalisation --Fig. 4.2 (a)

for t in target.unique():
    data_tmp = data[target==t]
    plt.scatter(data_tmp['petal length (cm)'], data_tmp['petal width (cm)'], label=t)
plt.legend()
plt.xlabel('petal length (cm)')
plt.ylabel('petal width (cm')
plt.xlim([1, 7])
plt.ylim([-1, 4])
plt.show()

Diagramme de dispersion après normalisation --Fig. 4.2 (b)

def standardize(X):
    return (X - np.mean(X)) / np.std(X)

std_data = standardize(data[['petal length (cm)', 'petal width (cm)']])
for t in target.unique():
    data_tmp = std_data[target==t]
    plt.scatter(data_tmp['petal length (cm)'], data_tmp['petal width (cm)'], label=t)
plt.legend()
plt.xlabel('petal length (cm)')
plt.ylabel('petal width (cm)')
plt.xlim([-2, 2])
plt.ylim([-2, 2])
plt.show()

Exemple d'exécution 4.2 Non corrélé

Rotation sur une matrice de vecteurs propres, c'est-à-dire décorrélation. Ci-dessous le code.

Diagramme de dispersion après décorrélation --Fig. 4.4 (b)

def no_correlation(X):
    cov_mat = np.cov(X, rowvar=0)
    eig_vals, eig_vecs = np.linalg.eig(cov_mat)
    return np.dot(X, eig_vecs)

no_corr_data = no_correlation(data[['petal length (cm)', 'petal width (cm)']])
for t in target.unique():
    data_tmp = no_corr_data[target==t]
    plt.scatter(data_tmp[:, 0], data_tmp[:, 1], label=t)
plt.legend()
plt.xlabel('petal length (cm)')
plt.ylabel('petal width (cm)')
plt.xlim([0, 8])
plt.ylim([-3, 3])
plt.show()

Dans cet exemple d'exécution, le signe de la composante y du vecteur propre est inversé par rapport à celui décrit dans le livre, de sorte que la direction de l'axe y du diagramme de dispersion est également inversée.

Sigma = np.cov(data[['petal length (cm)', 'petal width (cm)']], rowvar=0)
Lambda, S = np.linalg.eig(Sigma)
print('Lambda = \n', Lambda)
print('S = \n', S)
#production=>
# Lambda = 
#  [3.66123805 0.03604607]
# S = 
#  [[ 0.92177769 -0.38771882]
#  [ 0.38771882  0.92177769]]

Matrice de covariance après décorrélation

print('Sigma = \n', np.cov(no_corr_data, rowvar=0))
#production=>
# Sigma = 
#  [[3.66123805e+00 6.01365790e-16]
#  [6.01365790e-16 3.60460707e-02]]

Exemple d'exécution 4.3 Blanchiment

Le blanchiment est la centralisation, la décorrélation et la normalisation de l'écart type à 1. Ci-dessous le code.

Diagramme de diffusion après blanchiment --Fig.4.6 (b)

def whitening(X):
    cov_mat = np.cov(X, rowvar=0)
    eig_vals, eig_vecs = np.linalg.eig(cov_mat)
    diag_sqrt_eig_vals = np.diag(np.sqrt(eig_vals))

    return np.dot(np.dot(X-np.mean(X, axis=0), eig_vecs), np.linalg.inv(diag_sqrt_eig_vals.T))

whitened_data = whitening(data[['petal length (cm)', 'petal width (cm)']])
for t in target.unique():
    data_tmp = whitened_data[target==t]
    plt.scatter(data_tmp[:, 0], data_tmp[:, 1], label=t)
plt.legend()
plt.xlabel('petal length (cm)')
plt.ylabel('petal width (cm)')
plt.xlim([-3, 3])
plt.ylim([-3, 4])
plt.show()

Matrice de covariance après blanchiment

print('Sigma = \n', np.cov(whitened_data, rowvar=0))
#production=>
# Sigma = 
#  [[1.00000000e+00 1.64411249e-15]
#  [1.64411249e-15 1.00000000e+00]]

Exemple d'exécution 4.4 Fonction discriminante secondaire / fonction discriminante linéaire

Une fonction qui identifie une classe, en supposant que les probabilités conditionnelles de classe suivent une distribution normale. Ci-dessous le code.

--Lecture des données indiennes Pima

Les données sont empruntées à R datasets.

pima_train = pd.read_csv('https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/MASS/Pima.tr.csv')
pima_test = pd.read_csv('https://raw.githubusercontent.com/vincentarelbundock/Rdatasets/master/csv/MASS/Pima.te.csv')

Diagramme de dispersion des données d'entraînement --Fig.4.8

for t in ['Yes', 'No']:
    pima_tmp = pima_train[pima_train['type']==t]
    plt.scatter(pima_tmp['glu'], pima_tmp['bmi'], label=t)
plt.legend()
plt.xlabel('glu')
plt.ylabel('bmi')
plt.xlim([50, 200])
plt.ylim([15, 50])
plt.show()

--Mise en œuvre de la fonction d'identification secondaire

class QuadraticDiscriminantFunction:
    def __init__(self, threshold=0):
        self._S = None
        self._c = None
        self._F = None
        self._threshold = threshold
    
    def fit(self, X, Y):
        X_pos = X[Y==1]
        X_neg = X[Y==0]

        mu_pos = np.mean(X_pos, axis=0)
        mu_neg = np.mean(X_neg, axis=0)
        Sigma_pos = np.cov(X_pos, rowvar=0)
        Sigma_neg = np.cov(X_neg, rowvar=0)
        p_pos = len(Y[Y==1])/len(Y)
        p_neg = len(Y[Y==0])/len(Y)
        
        self._S = np.linalg.inv(Sigma_pos) - np.linalg.inv(Sigma_neg)
        self._c = np.dot(mu_neg, np.linalg.inv(Sigma_neg)) - np.dot(mu_pos, np.linalg.inv(Sigma_pos))
        self._F = np.dot(np.dot(mu_pos, np.linalg.inv(Sigma_pos)), mu_pos) - np.dot(np.dot(mu_neg, np.linalg.inv(Sigma_neg)), mu_neg)\
                    + np.log(np.linalg.det(Sigma_pos)/np.linalg.det(Sigma_neg))\
                    - 2*np.log(p_pos/p_neg)
    
    def predict(self, X):
        xSx = np.diag(np.dot(np.dot(X, self._S), X.T))
        cx = np.dot(self._c, X.T)
        Y_hat = xSx + 2*cx + self._F
        return np.where(Y_hat < self._threshold, 1, 0)

-Implémentation de la fonction discriminante linéaire

class LinearDiscriminantFunction:
    def __init__(self, threshold=0):
        self._c = None
        self._F = None
        self._threshold = threshold
    
    def fit(self, X, Y):
        X_pos = X[Y==1]
        X_neg = X[Y==0]

        mu_pos = np.mean(X_pos, axis=0)
        mu_neg = np.mean(X_neg, axis=0)
        Sigma_pos = np.cov(X_pos, rowvar=0)
        Sigma_neg = np.cov(X_neg, rowvar=0)
        p_pos = len(Y[Y==1])/len(Y)
        p_neg = len(Y[Y==0])/len(Y)
        Sigma_pool = p_pos*Sigma_pos + p_neg*Sigma_neg
        
        self._c = np.dot(mu_neg, np.linalg.inv(Sigma_pool)) - np.dot(mu_pos, np.linalg.inv(Sigma_pool))
        self._F = np.dot(np.dot(mu_pos, np.linalg.inv(Sigma_pool)), mu_pos) - np.dot(np.dot(mu_neg, np.linalg.inv(Sigma_pool)), mu_neg)\
                    - 2*np.log(p_pos/p_neg)
    
    def predict(self, X):
        cx = np.dot(self._c, X.T)
        Y_hat = 2*cx + self._F
        return np.where(Y_hat < self._threshold, 1, 0)

--Apprentissage et prédiction pour tracer des limites d'identification bayésiennes

X_train = pima_train[['glu', 'bmi']].values
Y_train = pima_train['type'].map({'Yes':1, 'No':0})

qdf = QuadraticDiscriminantFunction(threshold=0)
qdf.fit(X_train, Y_train)
ldf = LinearDiscriminantFunction(threshold=0)
ldf.fit(X_train, Y_train)

X1_min, X1_max = np.min(X_train[:, 0])-1, np.max(X_train[:, 0])+1
X2_min, X2_max = np.min(X_train[:, 1])-1, np.max(X_train[:, 1])+1
X1, X2 = np.meshgrid(np.arange(X1_min, X1_max, 1), np.arange(X2_min, X2_max, 0.5))
X = np.c_[X1.ravel(), X2.ravel()]

Y_qdf_pred = qdf.predict(X)
Y_ldf_pred = ldf.predict(X)

--Bordure de discrimination des baies obtenue à partir de la fonction discriminante secondaire --Fig. 4.9 (a)

plt.contour(X1, X2, Y_qdf_pred.reshape(X1.shape), colors='limegreen')
for t in ['Yes', 'No']:
    pima_tmp = pima_train[pima_train['type']==t]
    plt.scatter(pima_tmp['glu'], pima_tmp['bmi'], label=t)
plt.legend()
plt.xlabel('glu')
plt.ylabel('bmi')
plt.xlim([50, 200])
plt.ylim([15, 50])
plt.show()

Taux d'erreur de la fonction discriminante secondaire

X_test = pima_test[['glu', 'bmi']].values
Y_test = pima_test['type'].map({'Yes':1, 'No':0})

Y_train_qdf_pred = qdf.predict(X_train)
print('error rate (train) = {:.1f}%'.format(np.mean(Y_train_qdf_pred != Y_train)*100))
Y_test_qdf_pred = qdf.predict(X_test)
print('error rate (test) = {:.1f}%'.format(np.mean(Y_test_qdf_pred != Y_test)*100))

#production=>
# error rate (train) = 24.5%
# error rate (test) = 23.8%

--Bordure de discrimination des baies obtenue à partir de la fonction discriminante linéaire --Fig. 4.9 (b)

plt.contour(X1, X2, Y_ldf_pred.reshape(X1.shape), colors='limegreen')
for t in ['Yes', 'No']:
    pima_tmp = pima_train[pima_train['type']==t]
    plt.scatter(pima_tmp['glu'], pima_tmp['bmi'], label=t)
plt.legend()
plt.xlabel('glu')
plt.ylabel('bmi')
plt.xlim([50, 200])
plt.ylim([15, 50])
plt.show()

Taux d'erreur de la fonction discriminante linéaire

Y_train_ldf_pred = ldf.predict(X_train)
print('error rate (train) = {:.1f}%'.format(np.mean(Y_train_ldf_pred != Y_train)*100))
Y_test_ldf_pred = ldf.predict(X_test)
print('error rate (test) = {:.1f}%'.format(np.mean(Y_test_ldf_pred != Y_test)*100))
#production=>
# error rate (train) = 24.0%
# error rate (test) = 22.0%

--Courbe ROC de la fonction discriminante quadratique / fonction discriminante linéaire --Fig.4.10

qdf_tpr_list = []
qdf_fpr_list = []
ldf_tpr_list = []
ldf_fpr_list = []

for th in np.arange(-9, 7, 0.1):
    qdf = QuadraticDiscriminantFunction(threshold=th)
    qdf.fit(X_train, Y_train)
    Y_qdf_pred = qdf.predict(X_test)

    tpr_qdf = len(Y_test[(Y_test==0)&(Y_qdf_pred==0)]) / len(Y_test[Y_test==0])
    qdf_tpr_list.append(tpr_qdf)
    fpr_qdf = len(Y_test[(Y_test==1)&(Y_qdf_pred==0)]) / len(Y_test[Y_test==1])
    qdf_fpr_list.append(fpr_qdf)
                        
    ldf = LinearDiscriminantFunction(threshold=th)
    ldf.fit(X_train, Y_train)
    Y_ldf_pred = ldf.predict(X_test)
    
    tpr_ldf = len(Y_test[(Y_test==0)&(Y_ldf_pred==0)]) / len(Y_test[Y_test==0])
    ldf_tpr_list.append(tpr_ldf)
    fpr_ldf = len(Y_test[(Y_test==1)&(Y_ldf_pred==0)]) / len(Y_test[Y_test==1])
    ldf_fpr_list.append(fpr_ldf)

plt.plot(qdf_fpr_list, qdf_tpr_list, label='Quadratic')
plt.plot(ldf_fpr_list, ldf_tpr_list, label='Linear')
plt.legend()
plt.xlabel('false positive ratio')
plt.ylabel('true positive ratio')
plt.show()

en conclusion

J'ai essayé de reproduire l'exemple d'exécution du chapitre 4 de Hajipata en Python. Quant à l'exactitude du code, la figure dessinée correspond à peu près à la figure du livre, donc Yoshi! Nous vérifions uniquement le diplôme, donc si vous trouvez des erreurs, nous vous serions reconnaissants de bien vouloir faire une demande de modification.