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».
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'})
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]]
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()
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()
Rotation sur une matrice de vecteurs propres, c'est-à-dire décorrélation. Ci-dessous le code.
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]]
print('Sigma = \n', np.cov(no_corr_data, rowvar=0))
#production=>
# Sigma =
# [[3.66123805e+00 6.01365790e-16]
# [6.01365790e-16 3.60460707e-02]]
Le blanchiment est la centralisation, la décorrélation et la normalisation de l'écart type à 1. Ci-dessous le code.
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()
print('Sigma = \n', np.cov(whitened_data, rowvar=0))
#production=>
# Sigma =
# [[1.00000000e+00 1.64411249e-15]
# [1.64411249e-15 1.00000000e+00]]
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')
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()
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()
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()
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.
Recommended Posts