3.2 Classification et régression - Machine Learning en géosciences

Les problèmes qui appellent une réponse quantitative (une valeur numérique) relèvent de la régression ; ceux qui appellent une réponse qualitative (un booléen ou une catégorie) relèvent de la classification. Beaucoup de méthodes statistiques s’appliquent aux deux types de problèmes.

La classification binaire a deux classes de sortie. Elles finissent généralement en « A » et « non A ». Par exemple « séisme » ou « pas de séisme = bruit ». La classification multiclasse désigne celle qui compte plus de deux classes.

La classification suppose ici que nous connaissons les étiquettes : c’est une forme d’apprentissage supervisé.

1. Algorithmes de classification¶

Il existe plusieurs algorithmes de classification, que nous résumons ci-dessous avant de pratiquer. Pour chacun, posez trois questions dans cet ordre : à quelle question géoscientifique répond-il, quelle est l’idée statistique sous-jacente, et quand convient-il à vos données ? Les noms scikit-learn comptent le moins — un agent de code peut les fournir.

Régression logistique : cet échantillon d’eau est-il potable ? Cette forme d’onde est-elle un séisme ou du bruit ? La régression logistique trace une frontière droite dans l’espace des caractéristiques (features) et convertit la distance d’un échantillon à cette frontière en une probabilité entre 0 et 1. Parce qu’elle produit une probabilité et pas seulement une étiquette, elle convient aux problèmes où la réponse doit alimenter une prévision, et ses coefficients se lisent comme l’effet de chaque caractéristique. La leçon 3.6 lui est consacrée.
Analyse discriminante linéaire (LDA) : étant donné la densité et la susceptibilité magnétique, ce spécimen de roche est-il un granite ou un basalte ? La LDA traite chaque classe comme un nuage d’échantillons avec son propre centre mais une dispersion commune, et trouve la direction de l’espace des caractéristiques le long de laquelle les centres des classes se séparent le plus nettement, relativement à la dispersion interne de chaque classe ; la frontière de décision est une droite entre les nuages projetés. Là où la PCA (chapitre 2.12) cherche la direction de plus grande variance sans égard aux classes, la LDA cherche la direction qui distingue le mieux les classes. Parce qu’elle n’estime que les moyennes des classes et une dispersion commune, elle se comporte bien avec de petites collections étiquetées — quelques dizaines de spécimens par type de roche — là où les modèles flexibles surapprennent.
Bayésien naïf (Naive Bayes, NB) : lequel de trois types de roches a le plus plausiblement produit ces propriétés mesurées ? Le bayésien naïf demande, pour chaque classe, quelle serait la probabilité des valeurs de caractéristiques observées si l’échantillon appartenait à cette classe — en traitant chaque caractéristique comme si elle était indépendante des autres — et assigne la classe qui rend l’observation la plus probable. Cette hypothèse d’indépendance est rarement vraie des mesures géophysiques (température et humidité varient ensemble), et pourtant la méthode classe souvent bien, apprend sur très peu d’échantillons et n’a presque rien à régler.
K plus proches voisins (KNN) : si les échantillons d’entraînement les plus semblables à celui-ci sont surtout des basaltes, appelons-le basalte. KNN stocke les données d’entraînement et classe un nouvel échantillon par un vote parmi ses K voisins les plus proches dans l’espace des caractéristiques. « Le plus proche » dépend des unités — une densité en kg/m³ écraserait une susceptibilité mesurée en minuscules valeurs SI — de sorte que des caractéristiques en unités mélangées doivent d’abord être ramenées à des plages comparables. La méthode ne suppose rien de la forme de la frontière entre classes, mais elle faiblit quand le nombre de caractéristiques croît, et chaque prédiction exige une recherche dans l’ensemble d’entraînement stocké.
Machine à vecteurs de support (SVM) : ce sismogramme est-il un tir de carrière ou un séisme tectonique, quand il n’existe que quelques centaines d’exemples étiquetés à la main ? La SVM trouve la frontière qui sépare les classes avec le coussin le plus large — la marge — de sorte qu’un échantillon proche de la frontière devrait subir une grande perturbation pour changer de côté. Les noyaux (fonction de base radiale, polynomial) permettent à cette frontière de se courber. La SVM était la méthode de choix quand les catalogues géoscientifiques étiquetés comptaient des centaines d’événements plutôt que des millions, et elle reste forte dans ce régime des petits échantillons ; elle renvoie une décision, pas une probabilité.
Forêt aléatoire (RF) : cette combinaison de pente, de pluie et de lithologie marque-t-elle un versant comme propice aux glissements de terrain ? Une forêt aléatoire fait pousser de nombreux arbres de décision, chacun sur un sous-ensemble aléatoire des échantillons et des caractéristiques, et assigne la classe qui remporte le vote entre les arbres ; moyenner de nombreux arbres bruités réduit la variance de chacun d’eux. Elle tolère les caractéristiques en unités mélangées sans remise à l’échelle et indique quelles caractéristiques ont porté la classification — souvent la question géoscientifique qui compte. Ses fractions de vote ressemblent à des probabilités mais ne sont pas des probabilités calibrées — la leçon 3.6 y revient.
Réseaux de neurones artificiels (ANN) : quand la frontière entre les classes est trop emmêlée pour une droite, un coussin ou un vote — distinguer le trémor éruptif du bruit du vent dans un spectrogramme —, les réseaux de neurones composent des couches de fonctions simples en frontières de forme presque arbitraire. Cette flexibilité se paie en données : il leur faut beaucoup plus d’échantillons étiquetés qu’aux méthodes ci-dessus, et leurs décisions sont plus difficiles à interpréter. Le chapitre 4 les traite en profondeur.

Certains classifieurs gèrent nativement le multiclasse (descente de gradient stochastique — SGD ; classification par forêt aléatoire ; bayésien naïf). D’autres sont strictement binaires (régression logistique, machine à vecteurs de support — SVM). En pratique, scikit-learn enveloppe automatiquement les classifieurs binaires dans un schéma un-contre-tous (one-vs-rest), de sorte qu’ils fonctionnent malgré tout sur les problèmes multiclasses.

Exercice¶

Nous allons créer un jeu de données synthétique représentant trois types de roches : granite, basalte et grès. Chaque type aura des valeurs caractéristiques de densité et de susceptibilité magnétique.

🖥️ Diapositives du cours — Séance 12 (lun. 26 oct.)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.dummy import DummyClassifier
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
from sklearn.metrics import classification_report, confusion_matrix

Générons des données synthétiques


# Set random seed for reproducibility
np.random.seed(42)

# Number of samples per class
n_samples = 100

# Generate features for Granite
granite_density = np.random.normal(2.9, 0.2, n_samples)
granite_susceptibility = np.random.normal(0.0001, 0.0001, n_samples)
granite_label = ['Granite'] * n_samples

# Generate features for Basalt
basalt_density = np.random.normal(3.2, 0.2, n_samples)
basalt_susceptibility = np.random.normal(0.001, 0.0005, n_samples)
basalt_label = ['Basalt'] * n_samples

# Generate features for Sandstone
sandstone_density = np.random.normal(2.4, 0.2, n_samples)
sandstone_susceptibility = np.random.normal(0.00005, 0.00005, n_samples)
sandstone_label = ['Sandstone'] * n_samples

# Combine data
density = np.concatenate([granite_density, basalt_density, sandstone_density])
susceptibility = np.concatenate([granite_susceptibility, basalt_susceptibility, sandstone_susceptibility])
labels = np.concatenate([granite_label, basalt_label, sandstone_label])

# Create DataFrame
data = pd.DataFrame({
    'Density': density,
    'Magnetic Susceptibility': susceptibility,
    'Lithology': labels
})
data.head()

import seaborn as sns

sns.scatterplot(
    x='Density',
    y='Magnetic Susceptibility',
    hue='Lithology',
    data=data
)
plt.title('Rock Types Based on Density and Magnetic Susceptibility')
plt.show()

Découpons les données en trois : entraînement, validation et test (60/20/20). Chaque ensemble a un seul rôle. L’ensemble d’entraînement ajuste les paramètres du modèle. L’ensemble de validation compare les modèles et règle les hyperparamètres. L’ensemble de test n’est touché qu’une fois, tout à la fin, pour l’estimation finale et honnête de la performance.

X = data[['Density', 'Magnetic Susceptibility']]
y = data['Lithology']

# First split off the test set (20%), then split the rest into train (60%) and validation (20%)
X_temp, X_test, y_temp, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)
X_train, X_val, y_train, y_val = train_test_split(
    X_temp, y_temp, test_size=0.25, random_state=42
)
print(f"train: {len(X_train)}, validation: {len(X_val)}, test: {len(X_test)}")

train: 180, validation: 60, test: 60

baseline = DummyClassifier(strategy="most_frequent")
baseline.fit(X_train, y_train)
print(f"Majority-class baseline accuracy on validation set: {baseline.score(X_val, y_val):.3f}")

Majority-class baseline accuracy on validation set: 0.333

Entraînons un classifieur k-NN sur les données d’entraînement. Le nombre de voisins $k$ est un hyperparamètre : nous essayons plusieurs valeurs, notons chacune sur l’ensemble de validation et gardons la meilleure.

best_k, best_acc = None, 0.0
for k in [1, 3, 5, 9, 15]:
    clf = KNeighborsClassifier(n_neighbors=k)
    clf.fit(X_train, y_train)
    acc = clf.score(X_val, y_val)
    print(f"k = {k:2d}   validation accuracy = {acc:.3f}")
    if acc > best_acc:
        best_k, best_acc = k, acc
print(f"\nBest k on the validation set: {best_k}")

k =  1   validation accuracy = 0.683
k =  3   validation accuracy = 0.667
k =  5   validation accuracy = 0.700
k =  9   validation accuracy = 0.767
k = 15   validation accuracy = 0.750

Best k on the validation set: 9

Évaluons la performance sur l’ensemble de test, une seule fois, avec le $k$ retenu.

classifier = KNeighborsClassifier(n_neighbors=best_k)
classifier.fit(X_train, y_train)

y_pred = classifier.predict(X_test)
print(confusion_matrix(y_test, y_pred))
print(classification_report(y_test, y_pred))

[[12  4  0]
 [ 4 18  0]
 [ 1  3 18]]
              precision    recall  f1-score   support

      Basalt       0.71      0.75      0.73        16
     Granite       0.72      0.82      0.77        22
   Sandstone       1.00      0.82      0.90        22

    accuracy                           0.80        60
   macro avg       0.81      0.80      0.80        60
weighted avg       0.82      0.80      0.80        60

2. Algorithmes de régression¶

Les modèles de régression prédisent des valeurs numériques continues à partir des entrées.

Parmi les applications en géosciences : prédire les niveaux de nappe phréatique à partir d’enregistrements climatiques et de pompage, estimer des propriétés de sol à partir de mesures spectrales, et la descente d’échelle statistique de la température et des précipitations.

2.1 Régression linéaire¶

Soit $y$ les données et $\hat{y}$ la valeur prédite des données. Une régression linéaire générale se formule

$\hat{y} = w_0 + w_1 x_1 + ... + w_n x_n = h_w (\mathbf{x})$ .

$\mathbf{\hat{y}} = \mathbf{G} \mathbf{w}$ .

$y$ est un vecteur de données de longueur $m$ , $\mathbf{x}$ est un vecteur de caractéristiques de longueur $n$ . $\mathbf{w}$ est un vecteur de paramètres du modèle, et $h_w$ est appelée la fonction d’hypothèse, ou le modèle, de paramètres $w$ . Dans le cas le plus simple d’une régression linéaire en temps, la formulation devient :

$\hat{y} = w_0 + w_1 t$ ,

où $x_1 = t$ est la caractéristique temporelle.

Pour évaluer la performance du modèle, nous calculons un score de perte, ou résidu. Trois termes voisins ne recouvrent pas la même chose. Une fonction de perte mesure l’erreur sur un exemple. La fonction objectif est ce que l’optimiseur minimise : la perte agrégée sur les exemples, plus, le cas échéant, des termes de régularisation. Le mot coût est réservé ici aux conséquences réelles d’une erreur — une alerte manquée, une évacuation inutile — ou aux pondérations qu’on leur attribue. La perte la plus élémentaire est l’erreur quadratique moyenne (MSE, Mean Square Error) ; la régression linéaire de cet exercice n’ajoutant aucun terme de régularisation, la MSE moyennée sur l’ensemble d’entraînement est ici aussi la fonction objectif :

$MSE(\mathbf{x},h_w) = \frac{1}{m} \sum_{i=1}^{m} \left( h_w(\mathbf{x})_i - y_i \right)^2 = \frac{1}{m} \sum_{i=1}^{m} \left( \hat{y}_i - y_i \right)^2$ , dans le cas d’une régression linéaire.

L’équation normale est la solution de la régression linéaire qui minimise la MSE.

$\mathbf{w} = \left( \mathbf{x}^T\mathbf{x} \right)^{-1} \mathbf{x}^T \mathbf{y}$

À comparer avec le problème inverse classique posé par $\mathbf{d} = \mathbf{G} \mathbf{m}$ .

$\mathbf{m} = \left( \mathbf{G}^T\mathbf{G} \right)^{-1} \mathbf{G}^T \mathbf{d}$

On peut le résoudre avec le module d’algèbre linéaire de Numpy. Si $\left( \mathbf{x}^T\mathbf{x} \right)$ est singulière et ne peut pas être inversée, une matrice de rang inférieur appelée pseudo-inverse peut être calculée par décomposition en valeurs singulières. Nous avons aussi utilisé, lors d’une séance précédente, la fonction scikit-learn sklearn.linear_model.LinearRegression, qui est l’implémentation de la pseudo-inverse. Nous pratiquons ci-dessous ces inversions standard :

Autres algorithmes courants¶

Au-delà de la régression linéaire ordinaire, les algorithmes de régression courants incluent :

Régression polynomiale : à quelle vitesse le sol sèche-t-il quand la température monte, quand la réponse s’incurve au lieu de suivre une droite ? Ajouter les termes carrés et croisés des caractéristiques permet à la même machinerie des moindres carrés d’ajuster une courbe. Des degrés faibles suffisent en général ; les degrés élevés oscillent violemment hors de la plage des échantillons d’entraînement, ce qui est dangereux quand le modèle sera interrogé sur des conditions qu’il n’a jamais vues.
Régression à vecteurs de support (SVR) : la forme régressive de la SVM. Elle ajuste une fonction en ignorant les résidus plus petits qu’une tolérance choisie, si bien que seuls les échantillons situés hors de cette bande façonnent l’ajustement. Cela la rend indulgente au bruit de mesure et moins sensible à un capteur isolé qui part en pointe.
Régression par forêt aléatoire : quelle quantité de pluie est tombée, étant donné ces caractéristiques radar ? Un ensemble d’arbres de décision dont les prédictions sont moyennées ; il capture les relations courbes et les interactions entre caractéristiques, sans remise à l’échelle pour les unités mélangées. Une précaution : une forêt ne peut pas prédire au-delà de la plage de ses cibles d’entraînement — elle ne prévoira donc jamais une pluie plus forte que toutes celles sur lesquelles elle a été entraînée.
Arbres à gradient boosting (renforcement par gradient) : des arbres ajoutés séquentiellement, chacun ajustant les résidus laissés par les précédents. Sur les données tabulaires des géosciences — mesures de stations, diagraphies de forage, attributs de catalogues —, les arbres boostés sont souvent le choix le plus exact, au prix de plus d’hyperparamètres à régler et d’un plus grand appétit pour le surapprentissage qu’une forêt aléatoire.
Réseaux de neurones : ajustent des relations de forme presque arbitraire, avec assez d’échantillons. Pour les petits jeux de données structurés de ce chapitre, les méthodes plus simples ci-dessus les égalent en général ; le chapitre 4 précise quand ils prennent l’avantage.

Exercice¶

Objectif : prédire la teneur en eau du sol à partir de facteurs environnementaux.

Nous simulons une humidité du sol influencée par la température et l’humidité de l’air.

# Seeded random generator for reproducibility
rng = np.random.default_rng(42)

# Number of samples
n_samples = 300

# Generate environmental variables
temperature = rng.uniform(15, 35, n_samples)  # in degrees Celsius
humidity = rng.uniform(30, 90, n_samples)     # in percentage

# Generate soil moisture as a function of temperature and humidity,
# with a mild quadratic term in temperature
soil_moisture = (
    0.5 * humidity
    - 0.3 * temperature
    - 0.02 * (temperature - 25) ** 2
    + rng.normal(0, 2, n_samples)
)

# Create DataFrame
data = pd.DataFrame({
    'Temperature': temperature,
    'Humidity': humidity,
    'Soil Moisture': soil_moisture
})
data.head()

La fonction vraie comporte désormais une légère non-linéarité (le terme quadratique en température). C’est important pour l’exercice : la sélection de modèle n’a de sens que si la fonction vraie est inconnue. Si les données étaient exactement linéaires, la régression linéaire gagnerait par construction et comparer des modèles n’apprendrait rien.

Visualisation des données


fig = plt.figure(figsize=(10, 7))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(data['Temperature'], data['Humidity'], data['Soil Moisture'])
ax.set_xlabel('Temperature (°C)')
ax.set_ylabel('Humidity (%)')
ax.set_zlabel('Soil Moisture')
plt.tight_layout()
plt.show()

Découpons les données en ensembles d’entraînement, de validation et de test (60/20/20), avec les mêmes rôles que précédemment.

X = data[['Temperature', 'Humidity']]
y = data['Soil Moisture']

X_temp, X_test, y_temp, y_test = train_test_split(
    X, y, test_size=0.2, random_state=42
)
X_train, X_val, y_train, y_val = train_test_split(
    X_temp, y_temp, test_size=0.25, random_state=42
)
print(f"train: {len(X_train)}, validation: {len(X_val)}, test: {len(X_test)}")

train: 180, validation: 60, test: 60

Commençons par le modèle de référence trivial : prédire la moyenne des cibles d’entraînement pour chaque échantillon. Tout modèle de régression doit le battre.

from sklearn.dummy import DummyRegressor
from sklearn.metrics import mean_squared_error, r2_score

baseline = DummyRegressor(strategy="mean")
baseline.fit(X_train, y_train)
y_val_base = baseline.predict(X_val)
print(f"Baseline (mean prediction) on validation set: "
      f"MSE = {mean_squared_error(y_val, y_val_base):.2f}, "
      f"R2 = {r2_score(y_val, y_val_base):.3f}")

Baseline (mean prediction) on validation set: MSE = 87.15, R2 = -0.012

Puis le modèle réel le plus simple : la régression linéaire.

from sklearn.linear_model import LinearRegression

linreg = LinearRegression()
linreg.fit(X_train, y_train)
y_val_lin = linreg.predict(X_val)
mse_lin = mean_squared_error(y_val, y_val_lin)
print(f"Linear regression on validation set: "
      f"MSE = {mse_lin:.2f}, R2 = {r2_score(y_val, y_val_lin):.3f}")

Linear regression on validation set: MSE = 5.09, R2 = 0.941

Ensuite, un modèle polynomial de degré 2 : un pipeline scikit-learn (Pipeline) qui développe les caractéristiques en tous les termes de degré 2, puis ajuste une régression linéaire sur eux.

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures

polyreg = Pipeline([
    ("poly", PolynomialFeatures(degree=2)),
    ("linreg", LinearRegression()),
])
polyreg.fit(X_train, y_train)
y_val_poly = polyreg.predict(X_val)
mse_poly = mean_squared_error(y_val, y_val_poly)
print(f"Degree-2 polynomial on validation set: "
      f"MSE = {mse_poly:.2f}, R2 = {r2_score(y_val, y_val_poly):.3f}")

Degree-2 polynomial on validation set: MSE = 4.80, R2 = 0.944

Retenons le modèle dont la MSE de validation est la plus basse, puis rapportons sa performance une seule fois sur l’ensemble de test.

if mse_poly < mse_lin:
    best_name, best_model = "degree-2 polynomial", polyreg
else:
    best_name, best_model = "linear regression", linreg

y_test_pred = best_model.predict(X_test)
print(f"Selected model: {best_name}")
print(f"Test set: MSE = {mean_squared_error(y_test, y_test_pred):.2f}, "
      f"R2 = {r2_score(y_test, y_test_pred):.3f}")

Selected model: degree-2 polynomial
Test set: MSE = 4.02, R2 = 0.947

Le modèle polynomial gagne sur l’ensemble de validation parce que la fonction vraie contient un terme quadratique que le modèle linéaire ne peut pas représenter ; le score de test, utilisé une seule fois, confirme ce choix.