La classification binaire affecte chaque échantillon à l’une de deux classes. Notre tâche dans ce chapitre est un classique de la sismologie : étant donné une courte fenêtre de sismogramme, décider si elle contient un événement sismique (étiquette 1) ou seulement du bruit (étiquette 0). Plutôt que des formes d’onde brutes, nous travaillons avec quatre caractéristiques physiques calculées sur chaque fenêtre :
sta_lta: le rapport d’une amplitude moyenne à court terme à une amplitude moyenne à long terme. Une arrivée impulsive élève la moyenne à court terme et pousse le rapport bien au-dessus de 1. La distribution est à queue lourde et s’étend d’environ 0,4 à 40.kurtosis: le caractère impulsif de la distribution des amplitudes ; les arrivées en pointe l’augmentent.spectral_centroid_hz: la fréquence moyenne de la fenêtre pondérée par l’amplitude, environ 0,2 à 24 Hz. Les événements locaux portent plus d’énergie haute fréquence que le bruit microsismique d’origine océanique.dominant_freq_hz: la fréquence du pic spectral.
Pour les graphiques bidimensionnels de ce carnet, nous utilisons sta_lta et spectral_centroid_hz : une fenêtre d’événement tend à présenter à la fois un rapport STA/LTA élevé et un centroïde spectral élevé, tandis que le bruit se tient à de faibles valeurs des deux.
1. Données événement contre bruit¶
Nous chargeons un ensemble équilibré de 2 000 fenêtres, moitié événements et moitié bruit, depuis le paquet de données synthétiques du cours.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import mlgeo_synth
df = mlgeo_synth.detector_features(n=2000, event_fraction=0.5, seed=42)
df.head()fig, ax = plt.subplots()
for label, name, color in [(0, "noise", "tab:gray"), (1, "event", "tab:red")]:
sub = df[df["label"] == label]
ax.scatter(sub["sta_lta"], sub["spectral_centroid_hz"],
s=14, alpha=0.5, color=color, edgecolors="k", linewidths=0.3, label=name)
plt.xscale('log')
ax.set_xlabel("STA/LTA ratio")
ax.set_ylabel("Spectral centroid (Hz)")
ax.legend()
ax.grid(True)
plt.show()
L’axe des x en échelle logarithmique nous dit déjà comment construire la matrice de caractéristiques : prendre le logarithme d’un rapport à queue lourde comme STA/LTA est de l’ingénierie de caractéristiques standard, et nous le faisons avant de confier les caractéristiques à quelque modèle que ce soit.
X = np.column_stack([np.log10(df["sta_lta"]), df["spectral_centroid_hz"]])
y = df["label"].valuesAvant tout classifieur, le modèle de référence (baseline). Nous découpons d’abord les données : tout ce qui s’ajuste — une mise à l’échelle (scaler), un modèle de référence trivial, un classifieur — ne voit que l’ensemble d’entraînement, et l’ensemble de test reste intact jusqu’à l’évaluation. Le prédicteur trivial affecte chaque fenêtre à la classe la plus fréquente ; sur cet ensemble équilibré, cela donne une exactitude proche de 0,5, et chaque classifieur ci-dessous doit faire mieux.
from sklearn.dummy import DummyClassifier
from sklearn.model_selection import train_test_split
# split first: the baseline, like any model, is fit on the training split only
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.4, random_state=42)
baseline = DummyClassifier(strategy="most_frequent").fit(X_train, y_train)
print(f"Majority-class baseline accuracy: {baseline.score(X_test, y_test):.2f}")Majority-class baseline accuracy: 0.49
Nous commencerons par la méthode fondamentale : l’analyse discriminante linéaire (LDA).
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.preprocessing import StandardScaler
# define ML
clf = LinearDiscriminantAnalysis()
# normalize data: fit the scaler on the training split only, then apply it to
# both splits. Fitting on the full dataset would leak the test set's statistics
# into training — the exact flaw the audit exercise in 3.10 teaches you to catch.
scaler = StandardScaler().fit(X_train)
X_train = scaler.transform(X_train)
X_test = scaler.transform(X_test)
X_scaled = scaler.transform(X) # transform-only, for plotting the full set
# Fit the model.
clf.fit(X_train, y_train)
# calculate the mean accuracy on the given test data and labels.
score = clf.score(X_test, y_test)
print("The mean accuracy on the given test and labels is %f" % score)The mean accuracy on the given test and labels is 0.952500
from sklearn.inspection import DecisionBoundaryDisplay
ax = plt.subplot()
# plot the decision boundary as a background
DecisionBoundaryDisplay.from_estimator(clf, X_scaled, cmap='PiYG', alpha=0.8, ax=ax, eps=0.5)
ax.scatter(X_scaled[:, 0], X_scaled[:, 1], c=y, cmap='PiYG', alpha=0.6, edgecolors="k")
ax.set_xlabel("log10(STA/LTA), scaled")
ax.set_ylabel("Spectral centroid, scaled")
plt.show()
La LDA trace une seule droite dans l’espace des caractéristiques. Elle réussit bien ici parce que les deux classes sont à peu près séparables par une droite.
Essayons un autre classifieur : les K plus proches voisins (KNN).
from sklearn.neighbors import KNeighborsClassifier
# define ML
K = 5
clf = KNeighborsClassifier(K)
# same split, same train-fit scaling as above
# Fit the model.
clf.fit(X_train, y_train)
# calculate the mean accuracy on the given test data and labels.
score = clf.score(X_test, y_test)
print("The mean accuracy on the given test and labels is %f" % score)The mean accuracy on the given test and labels is 0.960000
# plot the decision boundary as a background
ax = plt.subplot()
DecisionBoundaryDisplay.from_estimator(clf, X_scaled, cmap='PiYG', alpha=0.8, ax=ax, eps=0.5)
ax.scatter(X_scaled[:, 0], X_scaled[:, 1], c=y, cmap='PiYG', alpha=0.6, edgecolors="k")
ax.set_xlabel("log10(STA/LTA), scaled")
ax.set_ylabel("Spectral centroid, scaled")
plt.show()
Des caractéristiques aux unités physiques différentes¶
Nous testons maintenant ce qui se passe quand nous ne normalisons pas les caractéristiques. Nos deux caractéristiques ont des sens physiques différents : le log-STA/LTA est sans dimension, le centroïde spectral est une fréquence. Supposons qu’un collègue nous transmette le centroïde en millihertz plutôt qu’en hertz. Rien de physique n’a changé, mais un axe est désormais numériquement mille fois plus grand.
# rebuild the raw features, spectral centroid now in millihertz
X_raw = np.column_stack([np.log10(df["sta_lta"]), df["spectral_centroid_hz"] * 1000.0])
# same KNN, but no scaling this time
clf_raw = KNeighborsClassifier(K)
Xr_train, Xr_test, yr_train, yr_test = train_test_split(X_raw, y, test_size=0.4, random_state=42)
clf_raw.fit(Xr_train, yr_train)
score = clf_raw.score(Xr_test, yr_test)
print("The mean accuracy without scaling is %f" % score)
# plot the decision boundary as a background
ax = plt.subplot()
DecisionBoundaryDisplay.from_estimator(clf_raw, X_raw, cmap='PiYG', alpha=0.8, ax=ax, eps=0.5)
ax.scatter(X_raw[:, 0], X_raw[:, 1], c=y, cmap='PiYG', alpha=0.6, edgecolors="k")
ax.set_xlabel("log10(STA/LTA)")
ax.set_ylabel("Spectral centroid (mHz)")
plt.show()The mean accuracy without scaling is 0.865000

L’exactitude chute et la frontière de décision devient un jeu de bandes horizontales : le classifieur n’utilise plus, en pratique, que le centroïde. Les classifieurs à base de distance ne voient que des nombres, pas des unités. La distance euclidienne dans le plan (log-STA/LTA, millihertz) est dominée par l’axe des millihertz, si bien que l’information STA/LTA est ignorée. La mise à l’échelle avec StandardScaler place les deux caractéristiques sur un pied d’égalité — c’est pourquoi nous avons normalisé avant l’ajustement.
2. Métriques de performance d’un classifieur¶
Dans un classifieur binaire, nous désignons l’une des deux classes comme positive, l’autre comme négative. Considérons N échantillons de données.
| Classe vraie \ Classe prédite | Négative | Positive | Total |
|---|---|---|---|
| Négative | Vrai négatif | Faux positif | n |
| Positive | Faux négatif | Vrai positif | p |
| Total | n’ | p’ | N |
Dans cet exemple, il y avait à l’origine étiquettes positives et étiquettes négatives. Nous obtenons au final échantillons prédits positifs et prédits négatifs.
Vrai positif TP : le nombre de données prédites positives qui étaient à l’origine positives.
Vrai négatif TN : le nombre de données prédites négatives qui étaient à l’origine négatives.
Faux positif FP : le nombre de données prédites positives mais qui étaient à l’origine négatives.
Faux négatif FN : le nombre de données prédites négatives mais qui étaient à l’origine positives.
Matrice de confusion :
On compte les cas où un élément de la classe A est classé dans la classe B : l’entrée est le nombre d’échantillons de classe vraie prédits en classe . scikit-learn trie les classes par étiquette, donc avec bruit = 0 (négatif) et événement = 1 (positif), la matrice qu’imprime confusion_matrix est :
La première ligne correspond aux fenêtres de bruit vrai, la seconde aux événements vrais ; les prédictions correctes se trouvent sur la diagonale. La matrice de confusion s’étend à la classification multi-classes : la matrice est alors KxK au lieu de 2x2. La meilleure matrice de confusion est proche d’une matrice diagonale, avec peu de termes hors diagonale.
Autres métriques de performance du modèle La performance du modèle peut être évaluée avec les quantités suivantes :
Erreur : la fraction des données mal classées
-> 0
Exactitude : la fraction des données correctement classées :
--> 1
Taux de vrais positifs : la fraction des échantillons réellement positifs qui sont correctement classés :
--> 1
Ce rapport est aussi le rappel (recall) ou la sensibilité.
Taux de vrais négatifs : la fraction des échantillons réellement négatifs qui sont correctement classés :
--> 1
Ce rapport est aussi la spécificité.
Précision : le rapport du nombre d’échantillons prédits dans la classe positive qui étaient effectivement positifs au nombre total d’échantillons prédits positifs.
--> 1
Score F1 :
--> 1.
La moyenne harmonique donne plus de poids au plus faible des deux, si bien que le score F1 n’est élevé que si le rappel et la précision le sont tous les deux.
Comment la précision et le rappel co-varient-ils ?
def returnPrecisionAndRecall(TP, FP, TN, FN):
precision = TP / (TP + FP) if (TP + FP) else 0
recall = TP / (TP + FN) if (TP + FN) else 0
return {'precision': precision, 'recall': recall}
# Okay, set how many true and false values are in the original dataset
actualTrueValues = 10
actualFalseValues = 10
sumValues = actualTrueValues + actualFalseValues
precisionRecallCollector = []
# Now, run a set of simulations
for n in range(5000):
TP = 0
FP = 0
TN = 0
FN = 0
# Begin by randomly setting the total number of true values returned
totalTrue = np.random.randint(0, high=sumValues)
# The total number of false values is sumValues - totalTrue
totalFalse = sumValues - totalTrue
# Partition totalTrue and totalFalse into TP, FP and TN, FN, respectively
if totalTrue != 0:
TP = np.random.randint(0, high=totalTrue)
FP = totalTrue - TP
if totalFalse != 0:
TN = np.random.randint(0, high=totalFalse)
FN = totalFalse - TN
thisPandR = returnPrecisionAndRecall(TP, FP, TN, FN)
precisionRecallCollector.append([thisPandR['precision'], thisPandR['recall']])pAndRArray = np.asarray(precisionRecallCollector)
plt.scatter(pAndRArray[:,0], pAndRArray[:,1], color='k', alpha=0.25)
plt.xlabel('Precision')
plt.ylabel('Recall')
Ce qui précède montre qu’on peut s’attendre à une large gamme de comportements entre précision et rappel. Ce n’est pas parce qu’un classifieur prédit bien une classe qu’il prédit bien l’autre. D’où l’importance de rapporter les deux métriques.
Affichons ces mesures pour notre classifieur KNN sur caractéristiques mises à l’échelle, évalué sur l’ensemble de test, avec scikit-learn.
from sklearn.metrics import confusion_matrix,precision_score,recall_score,f1_score
# Fit the model.
y_test_pred=clf.predict(X_test)
print("confusion matrix")
print(confusion_matrix(y_test,y_test_pred))
print("precison, recall")
print(precision_score(y_test,y_test_pred),recall_score(y_test,y_test_pred))
print("F1 score")
print(f1_score(y_test,y_test_pred))confusion matrix
[[384 12]
[ 20 384]]
precison, recall
0.9696969696969697 0.9504950495049505
F1 score
0.96
Un rapport complet et bien formaté de la performance s’obtient avec la fonction classification_report :
from sklearn.metrics import classification_report
print(f"Classification report for classifier {clf}:\n"
f"{classification_report(y_test, y_test_pred)}\n")Classification report for classifier KNeighborsClassifier():
precision recall f1-score support
0 0.95 0.97 0.96 396
1 0.97 0.95 0.96 404
accuracy 0.96 800
macro avg 0.96 0.96 0.96 800
weighted avg 0.96 0.96 0.96 800
Précision et rappel s’arbitrent l’un l’autre : augmenter la précision réduit le rappel.
Le classifieur utilise une valeur de seuil pour décider si une donnée appartient à une classe. Élever le seuil donne une meilleure précision, abaisser le seuil donne un meilleur rappel. Regardons les différentes valeurs de score.
Caractéristique opératoire du récepteur (ROC, Receiver Operating Characteristics)
Elle trace le taux de vrais positifs en fonction du taux de faux positifs. La courbe ROC est visuelle, mais on peut quantifier la performance du classifieur par l’aire sous la courbe (AUC). Idéalement, l’AUC vaut 1.

[source : https://
Nous calculons la courbe sur l’ensemble de test : les éditions antérieures de ce livre traçaient la ROC de l’ensemble d’entraînement, ce qui flatte le modèle.
y_scores = clf.predict_proba(X_test)
print(y_scores[:5])[[1. 0. ]
[1. 0. ]
[0. 1. ]
[0.8 0.2]
[0. 1. ]]
from sklearn.metrics import roc_curve
fpr,tpr,thresholds=roc_curve(y_test,y_scores[:,1])
plt.plot(fpr,tpr,linewidth=2);plt.grid(True)
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.plot([0,1],[0,1],'k--')
Nous explorons maintenant les différents classifieurs proposés par scikit-learn. Nous pouvons tester systématiquement leur performance et sauvegarder la précision, le rappel et le score F1.
Un mot de prudence sur une entrée : l’entraînement de GaussianProcessClassifier croît en avec le nombre d’échantillons. Cela passe ici à , mais c’est un piège sur les vrais catalogues de détection, qui comptent couramment 105 à 106 fenêtres.
from sklearn.svm import SVC
from sklearn.gaussian_process import GaussianProcessClassifier
from sklearn.gaussian_process.kernels import RBF
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier, AdaBoostClassifier
from sklearn.naive_bayes import GaussianNB
from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis
# define models
names = [
"Nearest Neighbors",
"Linear SVM",
"RBF SVM",
"Gaussian Process",
"Decision Tree",
"Random Forest",
"AdaBoost",
"Naive Bayes",
"QDA",
]
classifiers = [
KNeighborsClassifier(7),
SVC(kernel="linear", C=0.025),
SVC(gamma=2, C=1),
GaussianProcessClassifier(1.0 * RBF(1.0)),
DecisionTreeClassifier(max_depth=5),
RandomForestClassifier(max_depth=5, n_estimators=10, max_features=1),
AdaBoostClassifier(),
GaussianNB(),
QuadraticDiscriminantAnalysis(),
]3. Exploration des modèles¶
Explorez la performance de chacun de ces modèles sur les données synthétiques.
Sauvegardez dans un tableau les valeurs de précision, de rappel et de score F1.
Trouvez le modèle le plus performant
# reuse the train/test split and train-fit scaling from section 1
pre=np.zeros(len(classifiers))
rec=np.zeros(len(classifiers))
f1=np.zeros(len(classifiers))
for ii,iclass in enumerate(classifiers):
iclass.fit(X_train, y_train)
y_test_pred=iclass.predict(X_test)
pre[ii] =precision_score(y_test,y_test_pred)
rec[ii] =recall_score(y_test,y_test_pred)
f1[ii] =f1_score(y_test,y_test_pred)
df_scores=pd.DataFrame({'CLF name':names,'precision':pre,'recall':rec,'f1_score':f1})
print(df_scores) CLF name precision recall f1_score
0 Nearest Neighbors 0.969620 0.948020 0.958698
1 Linear SVM 0.969309 0.938119 0.953459
2 RBF SVM 0.967005 0.943069 0.954887
3 Gaussian Process 0.971939 0.943069 0.957286
4 Decision Tree 0.966837 0.938119 0.952261
5 Random Forest 0.962500 0.952970 0.957711
6 AdaBoost 0.959900 0.948020 0.953923
7 Naive Bayes 0.966921 0.940594 0.953576
8 QDA 0.966921 0.940594 0.953576
Détection déséquilibrée : quand l’exactitude ne veut plus rien dire¶
L’ensemble équilibré ci-dessus était une fiction de salle de classe bien commode. Sur une station réelle, les fenêtres d’événement sont rares : l’essentiel de la journée n’est que du bruit. Nous tirons maintenant 5 000 fenêtres avec une fraction d’événements de 2 pour cent, soit environ 100 événements pour 4 900 fenêtres de bruit — un déséquilibre de 1:50. Mêmes caractéristiques, même transformation logarithmique, même mise à l’échelle. Le découpage est stratifié, pour que la classe rare garde la même proportion dans l’entraînement et dans le test.
df_imb = mlgeo_synth.detector_features(n=5000, event_fraction=0.02, seed=42)
print(df_imb["label"].value_counts())
X_imb = np.column_stack([np.log10(df_imb["sta_lta"]), df_imb["spectral_centroid_hz"]])
y_imb = df_imb["label"].values
Xi_train, Xi_test, yi_train, yi_test = train_test_split(
X_imb, y_imb, stratify=y_imb, test_size=0.4, random_state=42
)
# same discipline as before: the scaler is fit on the training split only
scaler_imb = StandardScaler().fit(Xi_train)
Xi_train = scaler_imb.transform(Xi_train)
Xi_test = scaler_imb.transform(Xi_test)
baseline = DummyClassifier(strategy="most_frequent").fit(Xi_train, yi_train)
print(f"Majority-class baseline accuracy: {baseline.score(Xi_test, yi_test):.3f}")label
0 4900
1 100
Name: count, dtype: int64
Majority-class baseline accuracy: 0.980
Le prédicteur trivial qui ne déclare jamais d’événement atteint déjà 0,98. L’exactitude est désormais inutile comme métrique. Entraînons deux vrais classifieurs et regardons les métriques qui comptent pour la classe événement.
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import accuracy_score
logreg = LogisticRegression().fit(Xi_train, yi_train)
forest = RandomForestClassifier(n_estimators=100, random_state=42).fit(Xi_train, yi_train)
for name, model in [("Logistic regression", logreg), ("Random forest", forest)]:
yp = model.predict(Xi_test)
print(name)
print(f" accuracy : {accuracy_score(yi_test, yp):.3f}")
print(f" precision (event): {precision_score(yi_test, yp):.3f}")
print(f" recall (event) : {recall_score(yi_test, yp):.3f}")
print(f" F1 (event) : {f1_score(yi_test, yp):.3f}")Logistic regression
accuracy : 0.992
precision (event): 0.929
recall (event) : 0.650
F1 (event) : 0.765
Random forest
accuracy : 0.992
precision (event): 0.853
recall (event) : 0.725
F1 (event) : 0.784
Les deux modèles dépassent l’exactitude du modèle de référence de moins d’un point de pourcentage, et pourtant ils manquent environ un tiers des événements. Le rappel sur la classe événement est le nombre qui importe à l’analyste : il compte les séismes que le détecteur attrape réellement. Le compromis est opérationnel. Un détecteur réglé pour un rappel élevé inonde l’analyste de faux déclenchements à vérifier ; un détecteur réglé pour une précision élevée garde une liste de déclenchements propre mais manque les petits événements. Où se placer sur cette courbe est une décision qui porte sur le temps d’analyste et sur le coût scientifique, pas un détail de modélisation.
Nous pouvons voir tout le compromis d’un coup en traçant côte à côte la courbe ROC et la courbe précision-rappel.
from sklearn.metrics import PrecisionRecallDisplay, RocCurveDisplay
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(11, 4.5))
for name, model in [("Logistic regression", logreg), ("Random forest", forest)]:
RocCurveDisplay.from_estimator(model, Xi_test, yi_test, ax=ax1, name=name)
PrecisionRecallDisplay.from_estimator(model, Xi_test, yi_test, ax=ax2, name=name)
ax1.plot([0, 1], [0, 1], 'k--', linewidth=1)
ax1.set_title("ROC curve")
ax2.set_title("Precision-recall curve")
ax1.grid(True)
ax2.grid(True)
plt.tight_layout()
plt.show()
La courbe ROC a l’air excellente, et c’est exactement le problème. Avec un déséquilibre de 1:50, le taux de faux positifs divise par les quelque 1 960 fenêtres de bruit de l’ensemble de test, donc il reste minuscule même quand les faux positifs dépassent en nombre les vraies détections. La courbe précision-rappel divise au contraire par le nombre de positifs prédits, si bien qu’elle suit la charge de travail réelle de l’analyste : elle chute dès que la liste de déclenchements se remplit de bruit. Quand la classe positive est rare, lisez la courbe PR, pas la ROC.
Traiter le déséquilibre¶
Diagnostiquer le déséquilibre n’est que la moitié du travail. Deux leviers standard déplacent le point de fonctionnement sans collecter de nouvelles données :
class_weight="balanced"repondère la perte d’entraînement pour que chaque erreur sur un événement compte autant qu’environ 49 erreurs sur du bruit, ce qui pousse la frontière ajustée vers la classe rare.- Le déplacement de seuil :
predictcoupe la probabilité prédite à 0,5, mais rien n’impose ce choix. Abaisser le seuil échange de la précision contre du rappel le long de la courbe que nous venons de tracer.
logreg_bal = LogisticRegression(class_weight="balanced").fit(Xi_train, yi_train)
# threshold moving: same model as before, decision cut lowered from 0.5 to 0.25
proba = logreg.predict_proba(Xi_test)[:, 1]
yp_thresh = (proba >= 0.25).astype(int)
rows = []
for name, yp in [
("logreg, threshold 0.5", logreg.predict(Xi_test)),
("logreg, class_weight='balanced'", logreg_bal.predict(Xi_test)),
("logreg, threshold 0.25", yp_thresh),
]:
rows.append({"model": name,
"precision (event)": precision_score(yi_test, yp),
"recall (event)": recall_score(yi_test, yp),
"F1 (event)": f1_score(yi_test, yp)})
pd.DataFrame(rows).round(3)Les deux leviers élèvent le rappel sur les événements en payant en précision : la liste de déclenchements s’allonge, mais moins d’événements sont manqués. Aucun des deux n’est « meilleur » — ils glissent le long de la même courbe précision-rappel, et le choix du point de fonctionnement appartient à l’analyste, pas au modèle.
Conclusion¶
Le choix de la métrique est une décision de conception, pas un détail qu’on règle après coup. Il encode le coût opérationnel de chaque type d’erreur : un événement manqué coûte de la science, un faux déclenchement coûte du temps d’analyste. Décidez ce que coûte chaque erreur dans votre application, choisissez la métrique qui met un prix dessus, et alors seulement comparez les modèles.