Reproduire l'exemple d'exécution du chapitre 4 de Hajipata en Python

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'})
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()

4-2a.png

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()

4-2b.png

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.

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()

4-3b.png

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]]

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.

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()

4-6b.png

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')
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()

4-8.png

--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()

4-9a.png

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()

4-9b.png

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()

4-10.png

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.

Recommended Posts

Reproduire l'exemple d'exécution du chapitre 4 de Hajipata en Python
Reproduire l'exemple d'exécution du chapitre 5 de Hajipata en Python
Mesurons le résultat de l'exécution du programme avec C ++, Java, Python.
Vérifiez le comportement du destroyer en Python
[Python] PCA scratch dans l'exemple de "Introduction à la méthode d'analyse multivariée"
[Hikari-Python] Chapitre 07-02 Gestion des exceptions (exécution continue du programme par gestion des exceptions)
Principes de base pour exécuter NoxPlayer en Python
À la recherche du FizzBuzz le plus rapide en Python
Exemple pratique d'architecture hexagonale en Python
Un exemple de réponse à la question de référence de la session d'étude. Avec python.
Sortie du nombre de cœurs de processeur en Python
[Python] Trier la liste de pathlib.Path dans l'ordre naturel
Le contenu du didacticiel Python (chapitre 2) est résumé dans une puce.
Récupérer l'appelant d'une fonction en Python
Afficher le résultat du traitement de la géométrie en Python
Le contenu du didacticiel Python (chapitre 8) est résumé dans une puce.
Le contenu du didacticiel Python (chapitre 1) est résumé dans une puce.
Copiez la liste en Python
Le contenu du didacticiel Python (chapitre 10) est résumé dans une puce.
Découvrez la fraction de la valeur saisie en python
Trouvez la solution de l'équation d'ordre n avec python
L'histoire de la lecture des données HSPICE en Python
[Note] À propos du rôle du trait de soulignement "_" en Python
Résolution d'équations de mouvement en Python (odeint)
Sortie sous la forme d'un tableau python
Le contenu du didacticiel Python (chapitre 6) est résumé dans une puce.
Le contenu du didacticiel Python (chapitre 3) est résumé dans une puce.
le zen de Python
Comment passer le résultat de l'exécution d'une commande shell dans une liste en Python
Découvrez la bonne efficacité de calcul de la vectorisation en Python
Comment obtenir le nombre de chiffres en Python
[python] Récupère la liste des classes définies dans le module
Ruby, exécution de fragments de code Python de la sélection dans Emacs
L'histoire de FileNotFound en Python open () mode = 'w'
Apprenez le modèle de conception «Chaîne de responsabilité» en Python
Implémenter la solution de l'algèbre de Riccati en Python
Obtenir la taille (nombre d'éléments) de Union Find en Python
Ne pas être conscient du contenu des données en python
Utilisons les données ouvertes de "Mamebus" en Python
Implémentation de l'algorithme "Algorithm Picture Book" en Python3 (Heap Sort Edition)
[Python] Affiche toutes les combinaisons d'éléments de la liste
Obtenez l'URL de la destination de la redirection HTTP en Python
Un mémorandum sur la mise en œuvre des recommandations en Python
Pour faire l'équivalent de Ruby ObjectSpace._id2ref en Python
Vérifiez la nature atrophique de la distribution de probabilité en Python
[Python] Chapitre 01-02 À propos de Python (Exécution et installation de l'environnement de développement)
[Exemple d'amélioration de Python] Apprentissage des bases de Python sur un site gratuit en 2 semaines
[Apprentissage automatique] "Détection d'anomalies et détection de changement" Dessinons la figure du chapitre 1 en Python.
Vers la retraite de Python2
Trouver des erreurs en Python
Jugement d'équivalence d'objet en Python
Reproduire la méthode de division mutuelle euclidienne en Python
Exécution de commandes externes en Python
Implémentation du tri rapide en Python
À propos des fonctionnalités de Python
Le pouvoir des pandas: Python
Essayez de gratter les données COVID-19 Tokyo avec Python
Découvrez la largeur apparente d'une chaîne en python