Résumé: Exploration puis modélisation de données climatiques en utilisant Python et la librairie Scikit-learn. L'objectif est de prévoir pour le lendemain un possible dépassement d'un seuil de concentration en ozone à partir d'une prévision déterministe sur un maillage grossier et de variables climatiques locales. Estimation par différentes méthodes: régression logistique, k plus proches voisins, arbre de décision, agrégation de modèle, SVM. Comparaison des erreurs de prévision sur un échantillon test puis des courbes ROC. Itération sur plusieurs échantillons tests pour analyser la distribution de l'erreur de prévision. Ce calepin vient compléter l'étude faite avec R pour en comparer les deux approches.
Avertissement
L'objectif, sur ces données, est d'améliorer la prévision déterministe (MOCAGE), calculée par les services de MétéoFrance, de la concentration d'ozone dans certaines stations de prélèvement. Il s'agit d'un problème dit d'adaptation statistique ou post-traitement d'une prévision locale de modèles à trop grande échelle en s'aidant d'autre variables également gérées par MétéoFrance, mais à plus petite échelle (température, force du vent...).
La question posée reste: quelle est la meilleure stratégie pour prévoir l'occurrence d'un pic de pollution.
Comme avec R différentes méthodes sont testées : régression logistique, k plus proches voisins, arbre de décision, random forest, SVM. De façon générale on suppose que l'utilisateur dispose d'une installation python à jour. Le calepin a été testé avec la version 3.7.
Question subsidiaire quand préférer R ou Python ? Python conduit a des résultats (conclusions) identiques à ceux de R, moins complets pour leur interprétation, mais plus rapidement. Il s'agit des principales différences entre R pour "statisticien" et python pour "informaticien", on perd en interprétabilité mais on gagne en vitesse d'exécution.
Les données ont été extraites et mises en forme par le service concerné de Météo France. Elles sont décrites par les variables suivantes:
Ce sont des données "propres", sans trous, bien codées et de petites tailles. Elles présentent avant tout un caractère pédagogique.
Il est choisi ici de lire les données avec la librairie pandas pour bénéficier de la classe DataFrame. Ce n'est pas nécessaire pour l'objectif de prévision car les variables qualitatives ainsi construites ne peuvent être utilisées pour l'interprétation des modèles obtenus dans scikit-learn qui ne reconnaît pas la classe DataFrame.
In [ ]:
import pandas as pd
import numpy as np
# Lecture des données
## Charger les données ou les lire directement en précisant le chemin
path=""
# path="http://www.math.univ-toulouse.fr/~besse/Wikistat/data/"
ozone=pd.read_csv(path+"depSeuil.dat",sep=",",header=0)
# Vérification du contenu
ozone.head()
Ce qui suit permet d'affecter le bon type aux variables.
In [ ]:
ozone["STATION"]=pd.Categorical(ozone["STATION"],ordered=False)
ozone["JOUR"]=pd.Categorical(ozone["JOUR"],ordered=False)
ozone["O3obs"]=pd.DataFrame(ozone["O3obs"], dtype=float)
ozone.dtypes
In [ ]:
ozone.describe()
Même si les données ne présentent pas de défauts particuliers, une étude exploratoire préliminaire est indispensable afin de s'assurer le leur bonne cohérence, proposer d'éventuelles transformations et analyser les structures de corrélations ou plus généralement de liaisons entre les variables, de groupes des individus ou observations.
In [ ]:
%matplotlib inline
In [ ]:
import matplotlib.pyplot as plt
ozone["O3obs"].hist()
plt.show()
In [ ]:
ozone["MOCAGE"].hist()
plt.show()
Exercice Traiter ainsi toutes les variables. Ceci suggère des transformations pour une meilleure utilisation des modèles linéaires.
In [ ]:
from math import sqrt, log
ozone["SRMH2O"]=ozone["RMH2O"].map(lambda x: sqrt(x))
ozone["LNO2"]=ozone["NO2"].map(lambda x: log(x))
ozone["LNO"]=ozone["NO"].map(lambda x: log(x))
del ozone["RMH2O"]
del ozone["NO2"]
del ozone["NO"]
Exercice Vérifier l'opportunité de ces transformations (histogrammes des nouvelles variables).
Retirer les variables initiales et construire ci-dessous la variable "dépassement de seuil" pour obtenir le fichier qui sera effectivement utilisé.
In [ ]:
ozone["DepSeuil"]=ozone["O3obs"].map(lambda x: x > 150)
ozone.head()
In [ ]:
# scatter plot matrix des variables quantitatives
from pandas.plotting import scatter_matrix
scatter_matrix(ozone[["O3obs","MOCAGE","TEMPE","VentMOD","VentANG","SRMH2O","LNO2","LNO"]], alpha=0.2, figsize=(15, 15), diagonal='kde')
plt.show()
Q Commenter les relations entre les variables prises 2 à 2.
In [ ]:
from sklearn.decomposition import PCA
from sklearn.preprocessing import scale
# réduction des variables
X=scale(ozone[["O3obs","MOCAGE","TEMPE","VentMOD","VentANG","SRMH2O","LNO2","LNO"]])
Tous les résultats numétriques classiques sont fournis par l'implémentation de scikit-learn mais des efforts sont à produire pour construire les graphiques usuels généralement automatiquement produits par des librairies dédiées comme FactoMineR de R.
Les commandes suivantes permettent de réaliser une analyse en composantes principales sur les seules variables quantitatives. Par ailleurs la variable à modéliser (O3obs, concentration observée) n'est pas utilisée.
In [ ]:
pca = PCA()
## Estimation, calcul des composantes principales
C = pca.fit(X).transform(X)
## Décroissance de la variance expliquée
plt.plot(pca.explained_variance_ratio_)
plt.show()
In [ ]:
## distribution des composantes principales
plt.boxplot(C[:,0:20])
plt.show()
Q Commenter ces résultats: quel choix de la dimension?
Q Présence de valeurs atypiques.
In [ ]:
## Repésentation des individus
plt.figure(figsize=(5,5))
for i, j, nom in zip(C[:,0], C[:,1], ozone["DepSeuil"]):
color = "red" if nom else "blue"
plt.plot(i, j, "o",color=color)
plt.axis((-4,6,-4,6))
plt.show()
In [ ]:
## coordonnées et représentation des variables
coord1=pca.components_[0]*np.sqrt(pca.explained_variance_[0])
coord2=pca.components_[1]*np.sqrt(pca.explained_variance_[1])
fig = plt.figure(figsize=(5,5))
ax = fig.add_subplot(1, 1, 1)
for i, j, nom in zip(coord1,coord2, ozone[["O3obs","MOCAGE","TEMPE","VentMOD",
"VentANG","SRMH2O","LNO2","LNO"]].columns):
plt.text(i, j, nom)
plt.arrow(0,0,i,j,color='black')
plt.axis((-1.2,1.2,-1.2,1.2))
# cercle
c=plt.Circle((0,0), radius=1, color='gray', fill=False)
ax.add_patch(c)
plt.show()
Q Commenter la structure de corrélation des variables.
Q L'objectif est de définir une surface séparant les deux classes. Une discriminaiton linéaire (hyperplan) semble-t-elle possible?
Ce n'est pas utile ici mais une classification non supervisée est facile à obtenir à titre illustratif, par exemple en 4 classes, par l'algorithme k-means:
In [ ]:
from sklearn.cluster import KMeans
from sklearn.metrics import confusion_matrix
clust=KMeans(n_clusters=4)
clust.fit(X)
classe=clust.labels_
print(classe)
In [ ]:
## Repésentation des individus dans les coordonnées de l'acp.
plt.figure(figsize=(10,8))
plt.scatter(C[:,0], C[:,1], c=classe)
plt.show()
La recherche d'une meilleure méthode de prévision suit généralement le protocole suivant dont la première étape est déja réalisée.
Remarques
Transformation des données pour l'apprentissage.
Q Pourquoi les variables qualitatives sont-elles transformées en paquets d'indicatrices ou dummy variables?
Q Pourquoi le type data frame est transformé en une matrice.
In [ ]:
ozone.head()
In [ ]:
# Variables explicatives
ozoneDum=pd.get_dummies(ozone[["JOUR","STATION"]])
del ozoneDum["JOUR_0"]
ozoneQuant=ozone[["MOCAGE","TEMPE","VentMOD","VentANG","SRMH2O","LNO2","LNO"]]
dfC=pd.concat([ozoneDum,ozoneQuant],axis=1)
dfC.head()
In [ ]:
# variable à expliquer binaire
Yb=ozone["DepSeuil"].map(lambda x: int(x))
# variable à expliquer réelle
Yr=ozone["O3obs"]
In [ ]:
Yr.hist()
plt.show()
Extractions des échantillons d'apprentissage et test pour les deux types de modèles. Comme le générateur est initalisé de façon identique, ce sont les mêmes échantillons dans les deux cas.
In [ ]:
from sklearn.model_selection import train_test_split
X_train,X_test,Yb_train,Yb_test=train_test_split(dfC,Yb,test_size=200,random_state=11)
X_train,X_test,Yr_train,Yr_test=train_test_split(dfC,Yr,test_size=200,random_state=11)
L'étape suivante est une étape de standardisation des données ou normalisation. Les variables sont divisées par leur écart-type. Ce n'est pas utile dans le cas d'un modèle linéaire élémentaire car la solution est identique mais indispensbale pour beaucoup d'autres méthodes non linéaires (SVM, réseaux de neurones, modèles avec pénalisation). Cette étape est donc concrètement systématiquement exécutée pour éviter des soucis. Attention, les mêmes paramètres (moyennes, écarts-types) estimés sur l'échantillon d'apprentissage sont utilisés pour normaliser l'échantillon test.
In [ ]:
from sklearn.preprocessing import StandardScaler
# L'algorithme ds réseaux de neurones nécessite éventuellement une normalisation
# des variables explicatives avec les commandes ci-dessous
scaler = StandardScaler()
scaler.fit(X_train)
Xr_train = scaler.transform(X_train)
# Meme transformation sur le test
Xr_test = scaler.transform(X_test)
Les fonctions de modéles linéaires et linéaires généralisées sont limitées dans Scikit-learn et sans sorties numériques (tests) détaillées qui sont à rechercher dans une autre librairie (StatsModels). Dans les deux cas, les stratégies classiques (forward, backward, stepwise, Furnival et Wilson) de sélection de variables par optimisation d'un critère (Cp, AIC, BIC) ne semblent pas disponibles, même si AIC et BIC sont présents dans scikit-learn, et le type DataFrame (package pandas) n'est pas reconnu.
La façon efficace de procéder est donc d'introduire une pénalisation Lasso pour opérer une sélection de variables ou plutôt la sélection de variables quantitatives et d'indicatrices des modalités de celles qualitatives mais sans analyse fine des interactions comme cela est possible avec R.
Q Quel autre type de pénalisation est aussi utilisée en régression?
Q Quelle la méthode qui combine les deux?
A titre de comparaison, on trace la prévision de la concentration de l'échantillon test par la seule valeur du modèle Mocage ainsi que les résidus à ce modèle fonction de la valeur prédite (Mocage)
In [ ]:
plt.plot(X_train["MOCAGE"],Yr_train,"o")
plt.xlabel("Mocage")
plt.ylabel("O3 observee")
plt.show()
In [ ]:
from sklearn.metrics import r2_score
print("R2=",r2_score(Yr_train,X_train["MOCAGE"]))
In [ ]:
plt.plot(X_test["MOCAGE"],Yr_test,"o")
plt.xlabel("Mocage")
plt.ylabel("O3 observee")
plt.show()
In [ ]:
plt.plot(X_test["MOCAGE"],X_test["MOCAGE"]-Yr_test,"o")
plt.xlabel("Mocage")
plt.ylabel("Residus")
plt.show()
Q Commenter la qualité de ces résidus.
In [ ]:
# Erreur quadratique moyenne
from sklearn.metrics import mean_squared_error
print("MSE=",mean_squared_error(X_test["MOCAGE"],Yr_test))
In [ ]:
# Le coefficient de détermination
# peut être négatif en prévision avec un mauvais modèle,
# est nul si la prévision est constante égale à la moyennne
from sklearn.metrics import r2_score
print("R2=",r2_score(Yr_test,X_test["MOCAGE"]))
Comparer cette prévision déterministe (équation de Navier et Stockes) par l'adaptation statistique la plus élémentaire. Il s'agit d'une régression avec choix de modèle par régularisation avec une pénalisation lasso.
Q Quelles est la valeur par défaut du paramètre de pénalisation Lasso?.
In [ ]:
from sklearn import linear_model
regLasso = linear_model.Lasso()
regLasso.fit(Xr_train,Yr_train)
prev=regLasso.predict(Xr_test)
print("MSE=",mean_squared_error(Yr_test,prev))
In [ ]:
from sklearn.metrics import r2_score
print("R2=",r2_score(Yr_test,prev))
Le paramètre de pénalisation lasso est optimisé par validation croisée.
In [ ]:
from sklearn.model_selection import GridSearchCV
# grille de valeurs du paramètre alpha à optimiser
param=[{"alpha":[0.05,0.1,0.2,0.3,0.4,0.5,1]}]
regLasso = GridSearchCV(linear_model.Lasso(), param,cv=5,n_jobs=-1)
regLassOpt=regLasso.fit(Xr_train, Yr_train)
# paramètre optimal
regLassOpt.best_params_["alpha"]
print("Meilleur R2 = %f, Meilleur paramètre = %s" % (regLassOpt.best_score_,regLassOpt.best_params_))
Q Quelle validation croisée est exécutée?
Prévision avec la valeur optimale de alpha puis calcul et tracé des résidus.
In [ ]:
prev=regLassOpt.predict(Xr_test)
print("MSE=",mean_squared_error(prev,Yr_test))
print("R2=",r2_score(Yr_test,prev))
In [ ]:
plt.plot(prev,Yr_test,"o")
plt.xlabel(u"O3 Prédite")
plt.ylabel("O3 observee")
plt.show()
In [ ]:
plt.plot(prev,Yr_test-prev,"o")
plt.xlabel(u"Prédites")
plt.ylabel(u"Résidus")
plt.hlines(0,40,220)
plt.show()
Q Comparer ces résidus avec ceux précédents (mocage) et noter l'amélioration.
Q Commenter la forme du nuage et donc la validité du modèle.
L'interprétation nécessite de connaître les valeurs des coefficients du modèle alors que l'objet regLassOpt issu de GridSearchCV ne retient pas les paramètres estimés. Il faut donc le ré-estimer avec la valeur optimale du paramètre de pénalisation si l'on souhaite afficher ces coefficients.
In [ ]:
# Coefficients
regLasso=linear_model.Lasso(alpha=regLassOpt.best_params_['alpha'])
model_lasso=regLasso.fit(Xr_train,Yr_train)
model_lasso.coef_
In [ ]:
coef = pd.Series(model_lasso.coef_, index = X_train.columns)
print("Lasso conserve " + str(sum(coef != 0)) +
" variables et en supprime " + str(sum(coef == 0)))
In [ ]:
imp_coef = coef.sort_values()
plt.rcParams['figure.figsize'] = (8.0, 10.0)
imp_coef.plot(kind = "barh")
plt.title(u"Coefficients du modèle lasso")
Q Noter les conséquences de la pénalisation; interpréter l'effet de chaque variable sur la concentration en ozone.
C'est ici qu'apparaît une insuffisance de la librairie python. Il faudrait construire "à la main" ou utiliser la librairie Statsmodels pour afficher les statistiques des tests et p-valeurs. Même avec ces compléments, la prise en compte des interactions et de leur sélection ne sont pas prévues. De plus l'interprétation est compliquée par l'éclatement de chaque variable qualitative en paquets d'indicatrices. C'est encore compréhensible avec peu de variables mais devient rapidement inexploitable.
Le graphe quivant permet d'identifier les bonnes et mauvaises prévisions de dépassement du seuil légal, ici fixé à $ 150 \mu g $.
In [ ]:
plt.plot(prev,Yr_test,"o")
plt.xlabel(u"Valeurs prédites")
plt.ylabel(u"O3 observée")
plt.hlines(150,50,300)
plt.vlines(150,0,300)
plt.show()
In [ ]:
# Dénombrement des erreurs par
# matrice de confusion
table=pd.crosstab(prev>150,Yr_test>150)
print(table)
Q Observer l'asymétrie de cette matrice. A quoi est-elle due au moins en partie ?
Scikit-learn propose d'autres procédures d'optimisation du paramètre de régularisation lasso par validation croisée en régression; lassoCV utilise un algorithme de coordinate descent, sans calcul de dérivée puisque la norme l1 n'est pas dérivable, tandis que lassoLarsCV est basée sur l'algorithme de least angle regression. Ces fonctions permettent de tracer également les chemins de régularisation. Voici l'exemple de lassoCV qui offre plus d'options.
In [ ]:
from sklearn.linear_model import LassoCV, LassoLarsCV
model = LassoCV(cv=5, alphas=np.array(range(1,50,1))/20.,n_jobs=-1,random_state=13).fit(Xr_train,Yr_train)
m_log_alphas = -np.log10(model.alphas_)
plt.figure()
# ymin, ymax = 2300, 3800
plt.plot(m_log_alphas, model.mse_path_, ':')
plt.plot(m_log_alphas, model.mse_path_.mean(axis=-1), 'k',
label='MSE moyen', linewidth=2)
plt.axvline(-np.log10(model.alpha_), linestyle='--', color='k',
label='alpha: optimal par VC')
plt.legend()
plt.xlabel('-log(alpha)')
plt.ylabel('MSE')
plt.title('MSE de chaque validation: coordinate descent ')
plt.show()
Q Vérifier que c'est bien la même valeur optimale que celle précédemment trouvée.
Tracés des chemins de régularisation.
In [ ]:
from itertools import cycle
from sklearn.linear_model import lasso_path
alphas_lasso, coefs_lasso, _ = lasso_path(Xr_train,Yr_train, alphas=np.array(range(1,50,1))/20.,)
plt.figure()
ax = plt.gca()
styles = cycle(['-', '--', '-.', ':'])
neg_log_alphas_lasso = -np.log10(alphas_lasso)
for coef_l, s in zip(coefs_lasso, styles):
l1 = plt.plot(neg_log_alphas_lasso, coef_l, linestyle=s,c='b')
plt.xlabel('-Log(alpha)')
plt.ylabel('Coefficients')
plt.show()
La même démarche est déroulée mais en modélisant directement la variable binaire Yb de dépassement ou non du seuil. Il s'agit d'une régression logistique avec toujours une pénalisation Lasso pour opérer une sélection de variables.
In [ ]:
from sklearn.linear_model import LogisticRegression
In [ ]:
# Optimisation du paramètre de pénalisation
# grille de valeurs
param=[{"C":[1,1.2,1.5,1.7,2,3,4]}]
logit = GridSearchCV(LogisticRegression(penalty="l1",solver="liblinear"), param,cv=5,n_jobs=-1)
logitOpt=logit.fit(Xr_train, Yb_train) # GridSearchCV est lui même un estimateur
# paramètre optimal
logitOpt.best_params_["C"]
print("Meilleur score = %f, Meilleur paramètre = %s" % (1.-logitOpt.best_score_,logitOpt.best_params_))
In [ ]:
# erreur sur l'échantillon test
1-logitOpt.score(Xr_test, Yb_test)
Le modèle "optimal" obtenu est utilisé pour prédire l'échantillon test et estimer ainsi, sans biais, une erreur de prévision.
La matrice de confusion croise les dépassements de seuils prédits avec ceux effectivement observés.
In [ ]:
# Prévision
y_chap = logitOpt.predict(Xr_test)
# matrice de confusion
table=pd.crosstab(y_chap,Yb_test)
print(table)
L'interprétation du modèle est basée sur les valeurs des coefficients avec les mêmes difficultés ou restrictions que pour la régression. Attention, GridSearch ne retient pas les coefficients, il faut les ré-estimer.
In [ ]:
# Coefficients
logitLasso=LogisticRegression(penalty="l1",C=logitOpt.best_params_['C'],
solver="liblinear")
logitCoef=logitLasso.fit(Xr_train,Yb_train).coef_
print(logitCoef[0])
In [ ]:
coef = pd.Series(logitCoef[0], index = X_train.columns)
print("Lasso conserve " + str(sum(coef != 0)) +
" variables et en supprime " + str(sum(coef == 0)))
In [ ]:
imp_coef = coef.sort_values()
plt.rcParams['figure.figsize'] = (6.0, 6.0)
imp_coef.plot(kind = "barh")
plt.title(u"Coefficients du modèle lasso")
Q Interpréter l'effet des variables retenues.
In [ ]:
from sklearn.metrics import roc_curve
probas_ = LogisticRegression(penalty="l1", solver="liblinear",
C=logitOpt.best_params_['C']).fit(X_train, Yb_train).predict_proba(X_test)
fpr, tpr, thresholds = roc_curve(Yb_test, probas_[:,1])
plt.plot(fpr, tpr, lw=1)
plt.xlabel('Taux de faux positifs')
plt.ylabel('Taux de vrais positifs')
plt.show()
Q Commenter la courbe ROC à propos du choix de la valeur seuil.
Voici un cas d'application d'analyses discriminantes non paramétriques, celles paramétriques (gaussienes) linéaires et quadratiques sont également présentes dans scikit-learn mais laissées en exercice.
Le paramètre de compléxité ($k$) est optimisé sur une grille prédéfinie en minimisant l'erreur estimée par validation croisée; scikit-learn propose de nombreuses options de validation croisée.
In [ ]:
from sklearn.neighbors import KNeighborsClassifier
# Optimisation de k
# grille de valeurs
param_grid=[{"n_neighbors":list(range(1,15))}]
knn=GridSearchCV(KNeighborsClassifier(),param_grid,cv=5,n_jobs=-1)
knnOpt=knn.fit(Xr_train, Yb_train) # GridSearchCV est lui même un estimateur
# paramètre optimal
knnOpt.best_params_["n_neighbors"]
print("Meilleur score = %f, Meilleur paramètre = %s" % (1.-knnOpt.best_score_,knnOpt.best_params_))
In [ ]:
# Estimation de l'erreur de prévision sur l'échantillon test
1-knnOpt.score(Xr_test,Yb_test)
In [ ]:
# Prévision de l'échantillon test
y_chap = knnOpt.predict(Xr_test)
# matrice de confusion
table=pd.crosstab(y_chap,Yb_test)
print(table)
Exercice Compléter les résultats en utilisant la fonction KNeighborsRegressor pour modéliser la concentration; optimiser $k$, calculer la prévision de l'échantillon test, tracer le graphe des résidus, calculer le MSE sur l'échantillon test.
Les arbres binaires de décision : discrimination ou régression, sont bien implémentés dans scikit-learn mais avec une insuffisance pour leur élagage. Ce n'est pas une pénalisation de la complexité, et donc précisément le nombre de feuilles qui est optimisé, mais la profondeur globale de l'arbre au risque d'élaguer, à une profondeur donnée, des feuilles importantes ou de conserver des feuilles ambigües.
Comme précédemment, la validation croisée permet d'optimiser le paramètre sur une grille.
In [ ]:
from sklearn.tree import DecisionTreeClassifier
# Optimisation de la profondeur de l'arbre
param=[{"max_depth":list(range(2,10))}]
tree= GridSearchCV(DecisionTreeClassifier(),param,cv=10,n_jobs=-1)
treeOpt=tree.fit(Xr_train, Yb_train)
# paramètre optimal
print("Meilleur score = %f, Meilleur paramètre = %s" % (1. - treeOpt.best_score_,treeOpt.best_params_))
In [ ]:
# Estimation de l'erreur de prévision
1-treeOpt.score(Xr_test,Yb_test)
In [ ]:
# prévision de l'échantillon test
y_chap = treeOpt.predict(Xr_test)
# matrice de confusion
table=pd.crosstab(y_chap,Yb_test)
print(table)
Autre difficulté dans la représentation d'un arbre de décision binaire. Le logiciel conseillé (Graphviz) semble délicat d'installation et d'utilisation pour un néophyte. Il est possible de lister la construction des noeuds avec quelques lignes de commande.
In [ ]:
from sklearn.tree import export_graphviz
from sklearn.externals.six import StringIO
import pydotplus
treeG=DecisionTreeClassifier(max_depth=treeOpt.best_params_['max_depth'])
treeG.fit(Xr_train,Yb_train)
dot_data = StringIO()
export_graphviz(treeG, out_file=dot_data)
graph=pydotplus.graph_from_dot_data(dot_data.getvalue())
graph.write_png("treeOpt.png")
In [ ]:
from IPython.display import Image
Image(filename='treeOpt.png')
Q Que dire de l'interprétation de l'arbre? Comparer les rôles des variables avec le modèle logit.
Exercice Compléter les résultats en utilisant la fonction DecisionTreeRegressor pour modéliser concentration; optimiser la profondeur, calculer la prévision de l'échantillon test, tracer les résidus, calculer le MSE sur l'échantillon test.
Les réseaux neuronaux (perceptron multicouche) ne sont présents dans le package Scikit-learn qu'à partir de la version 0.18. Les méthodes profondes (deep learning) nécessitent l'installation des librairies theano et Lasagne ou theano, TensorFlow et Keras. Ces dernières sont nettement plus complexes à installer, surtout sous Windows. Elles feront l'objet d'un autre tutoriel.
In [ ]:
from sklearn.neural_network import MLPClassifier
Définition des paramètres dont le nombre de neurones et alpha qui règle la régularisation par défaut 10-5. Le nombre de neurones est optimisé mais ce peut être alpha avec un nombre grand de neurones. Le nombre max d'itérations par défaut (200) semble insuffisant. Il est fixé à 500.
In [ ]:
param_grid=[{"hidden_layer_sizes":list([(5,),(6,),(7,),(8,)])}]
nnet= GridSearchCV(MLPClassifier(max_iter=500),param_grid,cv=10,n_jobs=-1)
nnetOpt=nnet.fit(Xr_train, Yb_train)
# paramètre optimal
print("Meilleur score = %f, Meilleur paramètre = %s" % (1. - nnetOpt.best_score_,nnetOpt.best_params_))
In [ ]:
# Estimation de l'erreur de prévision sur le test
1-nnetOpt.score(Xr_test,Yb_test)
In [ ]:
# prévision de l'échantillon test
y_chap = nnetOpt.predict(Xr_test)
# matrice de confusion
table=pd.crosstab(y_chap,Yb_test)
print(table)
Exercice Remplacer ensuite la fonction MLPClassifier par celle MLPRegressor de régression. Optimiser le paramètre, calculer la prévision, les résidus, le MSE.
La librairie randomForest de R utilise le programme historique développé par Breiman et Cutler(2001) et interfacé par Liaw et Wiener. Cette interface est toujours mise à jour mais il n'est pas sûr que le programme original continue d'évoluer depuis 2004. Pour des tailles importantes d'échantillons, quelques milliers, cette implémentation atteint des temps d'exécution rédhibitoires (cf. cet exemple) au contraire de celle en Python dont gestion mémoire et capacité de parallélisation ont été finement optimisées par Louppe et al.(2014).
De même que le boosting, deux fonctions de forêt sont proposés dans scikit-learn ; une pour la régression et une pour la classification ainsi qu'une version "plus aléatoire". Par rapport à la version originale de R, moins d'options sont proposées mais l'utilisation de base est très similaire avec le même jeu de paramètres.
Q20 Identifier les paramètres, les valeurs par défaut.
In [ ]:
from sklearn.ensemble import RandomForestClassifier
# définition des paramètres
forest = RandomForestClassifier(n_estimators=500,
criterion='gini', max_depth=None,
min_samples_split=2, min_samples_leaf=1,
max_features='auto', max_leaf_nodes=None,
bootstrap=True, oob_score=True)
# apprentissage
rfFit = forest.fit(Xr_train,Yb_train)
print(1-rfFit.oob_score_)
Comparer l'erreur out-of-bag ci-dessus avec celle sur l'échantillon test.
In [ ]:
# erreur de prévision sur le test
1-rfFit.score(Xr_test,Yb_test)
Optimisation par validation croisée du nombre de variables tirés aléatoirement lors de la construction de chaque noeud.
In [ ]:
param=[{"max_features":list(range(2,10,1))}]
rf= GridSearchCV(RandomForestClassifier(n_estimators=100),
param,cv=5,n_jobs=-1)
rfOpt=rf.fit(Xr_train, Yb_train)
# paramètre optimal
print("Meilleur score = %f, Meilleur paramètre = %s" % (1. - rfOpt.best_score_,rfOpt.best_params_))
Plusieurs exécutions, rendues aléatoires par la validation croisée, peuvent conduire à des valeurs "optimales" différentes de ce paramètre sans pour autant nuire à la qualité de prévision sur l'échantillon test.
In [ ]:
# erreur de prévision sur le test
1-rfOpt.score(Xr_test,Yb_test)
Exercice Tester différentes valeurs de min_samples_split de celle trouvée optimale. Conclusion sur la sensibilité de l'optimisation de ce paramètre ?
In [ ]:
# prévision
y_chap = rfFit.predict(Xr_test)
# matrice de confusion
table=pd.crosstab(y_chap,Yb_test)
print(table)
Comme avec R, il est possible de calculer un indicateur d'importance des variables pour aider à une forme d'interprétation. Celui-ci dépend de la position de la variable dans l'arbre et correspond donc au mean decrease in Gini index de R plutôt qu'au mean descrease in accuracy. La forêt doit être réestimée car GridSearch ne connaît pas le paramètre d'importance.
In [ ]:
rf= RandomForestClassifier(n_estimators=100,max_features=2)
rfFit=rf.fit(Xr_train, Yb_train)
# Importance décroissante des variables
importances = rfFit.feature_importances_
indices = np.argsort(importances)[::-1]
for f in range(Xr_train.shape[1]):
print(dfC.columns[indices[f]], importances[indices[f]])
In [ ]:
# Graphe des importances
plt.figure()
plt.title("Importances des variables")
plt.bar(range(Xr_train.shape[1]), importances[indices])
plt.xticks(range(Xr_train.shape[1]), indices)
plt.xlim([-1, Xr_train.shape[1]])
plt.show()
Q Comparer les importances des variables et les sélections opérées précédemment.
Exercice Remplacer ensuite la fonction RandomForestClassifier par celle RandomForestRegressor de régression. Optimiser le paramètre, calculer la prévision, les résidus, le MSE.
Exercice Expérimenter également le boosting sur ces données en exécutant la fonction GradientBoostingClassifier opérant l'agorithme de gradient tree boosting.
Remarque: Une version "améliorée" de boosting mieux paralélisée et incluant d'autres paramètres (pénalisation), est proposé dans le package: XGBoost qui peut être utilisé à partir de Python mais aussi R, Julia ou Java. Nénamoins le choix est fait d'arrêter l'acharnement sur ces données; XGBoost est testé en python sur un autre jeu de données.
De nombreux paramètres sont associés à cette méthode. La liste est à consulter dans la documentation en ligne.
L'optimisation de la pénalisation (paramètre C) est recherchée sur une grille par validation croisée. Remarque: il serait nécessaire d'optimiser également la valeur du coefficient gamma lié au noyau gaussien ("écart-type").
Il est souvent nécessaire de normaliser des données avant d'opérer les SVM.
In [ ]:
from sklearn.svm import SVC
param=[{"C":[0.4,0.5,0.6,0.8,1,1.4]}]
svm= GridSearchCV(SVC(),param,cv=10,n_jobs=-1)
svmOpt=svm.fit(Xr_train, Yb_train)
# paramètre optimal
print("Meilleur score = %f, Meilleur paramètre = %s" % (1. - svmOpt.best_score_,svmOpt.best_params_))
In [ ]:
# erreur de prévision sur le test
1-svmOpt.score(Xr_test,Yb_test)
In [ ]:
# prévision de l'échantillon test
y_chap = svmOpt.predict(Xr_test)
# matrice de confusion
table=pd.crosstab(y_chap,Yb_test)
print(table)
Exercice Comme précédemment, remplacer ensuite la fonction SVC par celle SVR de régression. Optimiser le paramètre, calculer la prévision, les résidus; le MSE.
Dans toute méthode, la prévision de dépassement ou non est associée au choix d'un seuil qui est par défaut 0.5. L'optimisaiton de ce seuil dépend des coûts respectifs associés aux faux positifs et aux faux négatifs qui ne sont pas nécessairement égaux. La courbe ROC permet de représenter l'influence de ce seuil sur les taux de faux positifs et vrais positifs.
In [ ]:
from sklearn.metrics import roc_curve
listMethod=[["RF",rfOpt],["NN",nnetOpt],["Tree",treeOpt],["K-nn",knnOpt],["Logit",logitOpt]]
In [ ]:
for method in enumerate(listMethod):
probas_ = method[1][1].fit(Xr_train, Yb_train).predict_proba(Xr_test)
fpr, tpr, thresholds = roc_curve(Yb_test, probas_[:,1])
plt.plot(fpr, tpr, lw=1,label="%s"%method[1][0])
plt.xlabel('Taux de faux positifs')
plt.ylabel('Taux de vrais positifs')
plt.legend(loc="best")
plt.show()
Q22 Le critère d'AUC (aire sous la courbe) permet-il d'ordonner les courbes et donc les méthodes?
C'est à un taux de faux positif admissible et donc à valeur de seuil fixé qu'il faut choisir la méthode d'apprentissage à privilégier.
L'échantillon test est de taille modeste et donc l'estimation de l'erreur de prévision peut présenter une variance importante. Celle-ci est réduite en opérant une forme de validation croisée (Monte Carlo) en tirant plusieurs couples d'échantillon apprentissage et test pour itérer les traitements précédents. Les données sont normalisées pour toutes les méthodes car les autres que SVM et NN ne sont pas affectées.
Les fonctionnalités de scikit-learn se prètent bien à l'automatisation de ces traitements enchaînant extraction d'échantillons, estimation de plusieurs modèles, optimisation de leurs paramètres et estimation de l'erreur de prévision sur le test.
Le code est compact et d'exécution efficace car bien parallélisé par les fonctions utilisées.
In [ ]:
from sklearn.utils import check_random_state
import time
check_random_state(13)
tps0=time.perf_counter()
# définition des estimateurs
logit= LogisticRegression(penalty="l1",solver="liblinear")
knn = KNeighborsClassifier()
tree = DecisionTreeClassifier()
nnet = MLPClassifier(max_iter=600)
rf = RandomForestClassifier(n_estimators=100)
svm = SVC()
# Nombre d'itérations
B=3 # pour exécuter après le test, mettre plutôt B=30
# définition des grilles de paramètres
listMethGrid=[[svm,{"C":[0.4,0.5,0.6,0.8,1,1.4]}],
[rf,{"max_features":list(range(2,10,2))}],
[nnet,{"hidden_layer_sizes":list([(5,),(6,),(7,),(8,)])}],
[tree,{"max_depth":list(range(2,10))}],
[knn,{"n_neighbors":list(range(1,15))}],
[logit,{"C":[0.5,1,5,10,12,15,30]}]]
# Initialisation à 0 des erreurs pour chaque méthode (colonne) et chaque itération (ligne)
arrayErreur=np.empty((B,6))
for i in range(B): # itérations sur B échantillons test
# extraction apprentissage et test
X_train,X_test,Yb_train,Yb_test=train_test_split(dfC,Yb,test_size=200)
scaler = StandardScaler()
scaler.fit(X_train)
X_train = scaler.transform(X_train)
# Meme transformation sur le test
X_test = scaler.transform(X_test)
# optimisation de chaque méthode et calcul de l'erreur sur le test
for j,(method, grid_list) in enumerate(listMethGrid):
methodGrid=GridSearchCV(method,grid_list,cv=10,n_jobs=-1,
iid="TRUE").fit(X_train, Yb_train)
methodOpt = methodGrid.best_estimator_
methFit=methodOpt.fit(X_train, Yb_train)
arrayErreur[i,j]=1-methFit.score(X_test,Yb_test)
tps1=time.perf_counter()
print("Temps execution en mn :",(tps1 - tps0))
dataframeErreur=pd.DataFrame(arrayErreur,columns=["SVM","RF","NN","Tree","Knn","Logit"])
In [ ]:
# Distribution des erreurs de prévisions
# Les SVM présentant des erreurs atypiques sont laissés de côté.
dataframeErreur[["SVM","RF","NN","Tree","Knn","Logit"]].boxplot(return_type='dict')
plt.show()
In [ ]:
# Moyennes
dataframeErreur.mean()
Cet exemple, traité en R puis en Python, résume bien l'intérêt et le contexte des méthodes d'apprentissage.
Remarque Il est possible d'exécuter directement l'épisode 5 sans passer par toutes les étapes de classification supervisée. Il suffit d'exécuter jusqu'à la section 4.1 de l'épidode 1, phase exploratoire et préparation des échantillons, afin de construire les données utilisées dans les sections 12 et 13 d'imputation des données manquantes et de détection d'atypiques.
Les vraies données sont le plus souvent mitées par l'absence de données, conséquences d'erreurs de saisie, de pannes de capteurs... Les librairies de Python (pandas) offrent des choix rudimentaires pour faire des imputations de données manquantes quand celles-ci le sont de façon complètement aléatoire.
Le calepin R d'analyse de ces mêmes données propose une comparaison assez détaillée de deux stratégiées afin d'évaluer leurs performances respectives.
La première stratégie commence par imputer les données manquantes en les prévoyant par l'algorithme missForest. Une fois les données manquantes imputées, différentes méthodes de prévision sont utilisables comme précédemment. Deux sont exécutées: forêts aléatoires et extrem gradient boosting.
La deuxième stratégie évite l'étape d'imputation en exécutant directement un algorithme de prévision tolérant des données manquantes. Peu le fond, c'est le cas de XGBoost.
Sur ces données, mais sans gros effort d'optimisation de XGBoost, la première stratégie enchaînant missForest puis randomForest conduit à de meilleurs résultats. Seule celle-ci est employée dans ce calepin mais, bien évidemment, l'exécution de xgboost sans imputation préalable est une option également possible en Python.
Bien moins de méthodes sont proposées en Python, SCikit-learn ne proposant que des imputations basiques par la moyenne ou la médiane comme dans pandas. Néanmoins une imputation par prévision utilisant k-n, l'ACP ou des forêts aléatoires est disponible dans la librairie predictive_imputer.
Les commandes ci-dessous font appel aux fichiers suivants:
X données complètes initiales Xna les données avec des trous, XnaImp les données avec imputations ozoneLes données initiales de la base ozone sont reprises. Seule la variable à expliquer de dépassement de seuil est conservée. La première opération consiste à générer aléatoirement un certain taux de données manquantes par la fonction définie ci-dessous.
In [ ]:
import numpy as np
import numpy.ma as ma
import random
def input_nan(x, tx):
"""
x : a 2D matrix of float dtype
tx: the rate of nan value to put in the matrix
"""
n_total = x.shape[0] * x.shape[1]
mask = np.array([random.random() for _ in range(n_total)]).reshape(x.shape)<tx
mx = ma.masked_array(x, mask=mask, fill_value=np.nan)
return mx.filled()
In [ ]:
# données initiales avec
X=dfC
# Génération de 10% de valeurs manquantes
Xna=input_nan(X, .1)
In [ ]:
from predictive_imputer import predictive_imputer
# Choix du mode de prévision par *k*-pp, ACP ou forêt aléatoire
imputer = predictive_imputer.PredictiveImputer(f_model="RandomForest")
# Imputation
XnaImp = imputer.fit(Xna).transform(Xna.copy())
In [ ]:
# Données sans trous
X_train,X_test,Yb_train,Yb_test=train_test_split(X,Yb,test_size=200,random_state=11)
XnaImp_train,XnaImp_test,Yb_train,Yr_test=train_test_split(XnaImp,Yb,test_size=200,random_state=11)
In [ ]:
from sklearn.ensemble import RandomForestClassifier
# prévision sans trous
forest = RandomForestClassifier(n_estimators=500)
# apprentissage
rfFit = forest.fit(X_train,Yb_train)
# erreur de prévision
# erreur de prévision sur le test
1-rfFit.score(X_test,Yb_test)
In [ ]:
# prévision avec trous imputés
forest = RandomForestClassifier(n_estimators=500)
# apprentissage
rfFit = forest.fit(XnaImp_train,Yb_train)
# erreur de prévision
# erreur de prévision sur le test
1-rfFit.score(XnaImp_test,Yb_test)
Q Que dire de la qualité de prévision avec 10% de trous
Exercice Faire varier le taux de trous et étudier la dégradation de la prévision.
Exercice Comparer avec une approche directe de la prévision avec XGBoost sans imputation préalable.
Remarque Il est possible d'exécuter directement l'épisode 5 sans passer par toutes les étapes de classification supervisée. Il suffit d'exécuter jusqu'à la section 4.1 de l'épidode 1, phase exploratoire et préparation des échantillons, afin de construire les données utilisées dans les sections 6 et 7 d'imputation des données manquantes et de détection d'atypiques.
Le calepin R d'analyse de ces mêmes données propose une comparaison assez détaillée des scores de détection des anomalies. Comme dans R, Scikit-learn propose des fonctions en Pyhton de détection d'atypiques multidimensionnels. Les principales sont LOF et Isolation Forest dont les résultats sont comparés ci-dessous.
Les données sont restreintes aux seules variables quantitatives explicatives.
Q Quel est le rôle du paramètre k ci-dessous?
In [ ]:
ozone.head()
In [ ]:
ozoneR=ozone[["MOCAGE","TEMPE","VentMOD","VentANG","SRMH2O","LNO2","LNO"]]
ozoneR.head()
In [ ]:
from sklearn.neighbors import LocalOutlierFactor
clf = LocalOutlierFactor(n_neighbors=20, contamination=0.05 ) # choix de n_n par défaut
scoreLOF=clf.fit_predict(ozoneR)
scoreAtyp=-clf._decision_function(ozoneR)# opposé du LOF
plt.boxplot(scoreAtyp)
plt.show()
Q Comment se comporte le LOF en fonction de k?
Q Quel taux d'observations par défaut sont considérées comme atypiques?
In [ ]:
atypLofInd = clf.fit_predict(X)
L'analyse en composante principale est utilisée pour représenter les observations atypiques.
In [ ]:
## Repésentation des atypiques
plt.figure(figsize=(5,5))
for i, j, nom in zip(C[:,0], C[:,1], atypLofInd):
color = "red" if nom!=1 else "blue"
plt.plot(i, j, "o",color=color)
plt.axis((-4,6,-4,6))
plt.show()
In [ ]:
from sklearn.svm import OneClassSVM
clf=OneClassSVM(nu=0.1, gamma=0.01)
scoreSVM=clf.fit(ozoneR)
scoreAtypSVM=clf._decision_function(ozoneR)
plt.boxplot(scoreAtypSVM)
plt.show()
Q Quel taux d'atypiques par défaut?
In [ ]:
atypSVMInd = clf.predict(ozoneR)
In [ ]:
## Repésentation des atypiques
plt.figure(figsize=(5,5))
for i, j, nom in zip(C[:,0], C[:,1], atypSVMInd):
color = "red" if nom!=1 else "blue"
plt.plot(i, j, "o",color=color)
plt.axis((-4,6,-4,6))
plt.show()
In [ ]:
from sklearn.ensemble import IsolationForest
clf = IsolationForest(max_samples=1000, contamination=0.05,behaviour="new")
scoreIF=clf.fit(ozoneR)
scoreAtypIF=clf.decision_function(ozoneR)
plt.boxplot(scoreAtypIF)
plt.show()
In [ ]:
atypIFInd = clf.predict(ozoneR)
## Repésentation des atypiques
plt.figure(figsize=(5,5))
for i, j, nom in zip(C[:,0], C[:,1], atypIFInd):
color = "red" if nom!=1 else "blue"
plt.plot(i, j, "o",color=color)
plt.axis((-4,6,-4,6))
plt.show()
Q Les observations définies comme des anomalies se retrouve-t-elles généralement d'une approche à l'autre?
Remarques
randomForest de R, il ne semble pas exister de librairie proposant la détection d'anomalies relativemement à la construction d'un modèle de prévision y=f(X) par forêt aléatoire. Il importe de suivre l'évolution des librairies en cours de développement.
In [ ]:
In [ ]: