Un nouvel estimateur basé sur scikit-learn¶

On souhaite construire un estimateur qui estime une régression linéaire à coefficients positifs, une autre avec des coefficients uniquement négatifs puis pour finir une dernière régression linéaire qui considère les deux premières comme features.

Une régression linéaire minimise l’erreur $\sum_i \left\Vert X_i\theta - y_i \right\Vert^2$ . Le gradient est $\sum_i X_i'\left( X_i\theta - y_i \right)$ .

Comme le modèle souhaité est équivalent à une optimisation sous contrainte, on propose de le résoudre comme ceci :

On applique une itération de l’algorithme de la descente de gradient : $\theta_{t+1} = \theta_t - \epsilon_t \sum_i X_i'\left( X_i\theta - y_i \right)$ .
On ne garde que les coefficients positifs : $\theta_{t+1} = \max(0, \theta_t)$ .
On retourne à l’étape 1 ou on s’arrête si l’algorithme a convergé.

Le nouveau régresseur linéaire¶

[45]:

import sklearn
import matplotlib.pyplot as plt
import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.base import BaseEstimator, RegressorMixin, TransformerMixin
from sklearn.compose import ColumnTransformer, make_column_selector
from sklearn.ensemble import StackingRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from teachpyx.datasets.wines import load_wines_dataset

[14]:

class PositiveOrNegativeLinearRegression(RegressorMixin, BaseEstimator):
    """
    Apprends une régression linéaire avec des coefficients du même signe.

    :param epsilon: gradient step
    :param max_iter: number maximum of iterations
    :param positive: only positive weights (or negative if False)
    """

    def __init__(
        self, epsilon: float = 1.0, max_iter: int = 200, positive: bool = True
    ):
        super().__init__()
        self.epsilon = epsilon
        self.max_iter = max_iter
        self.positive = positive

    def fit(self, X, y):
        """
        Trains. L'algorithme d'apprentissage utilise le gradient global
        et restreint les coefficients à être positifs. On garde les derniers coefficients
        non tous nuls.
        """
        theta = np.random.randn(X.shape[1])
        theta = np.maximum(theta, 0) if self.positive else np.minimum(theta, 0)
        n_refresh = 0
        theta0 = None
        for i in range(self.max_iter):
            epsilon = self.epsilon * 10 / (i + 10)
            grad = X.T @ (X @ theta - y)
            theta = theta - epsilon * grad / X.shape[0]
            theta = np.maximum(theta, 0) if self.positive else np.minimum(theta, 0)
            if np.abs(theta).max() > 0:
                theta0 = theta
            if (
                np.abs(theta).max() == 0
                and i < self.max_iter - 5
                and n_refresh <= self.max_iter // 5
            ):
                # If all coefficients are null. The algorithm may be stuck.
                theta = np.random.randn(X.shape[1]) * epsilon
                theta = np.maximum(theta, 0) if self.positive else np.minimum(theta, 0)
                n_refresh += 1

        self.theta_ = theta if np.abs(theta).max() or theta0 is None else theta0
        return self

    def predict(self, X):
        "Predicts."
        if hasattr(X, "values"):
            # Input is a dataframe.
            return X.values @ self.theta_
        return X @ self.theta_

L’ordre d’héritage RegressorMixin, BaseEstimator est important pour que le modèle hérite des Tags. S’ils sont mal spécifiés, des supers modèles comme StackingRegressor créent une erreur. On vérifie que le modèle fonctionne.

[15]:

X = np.random.randn(5, 2)
y = X[:, 0] * 0.3 + X[:, 1] * 0.6 + np.random.randn(5) / 100
lr = LinearRegression().fit(X, y)
plr = PositiveOrNegativeLinearRegression(max_iter=50, positive=True)
plr.fit(X, y)
lr.coef_, plr.theta_

[15]:

(array([0.28899397, 0.59021441]), array([0.28446758, 0.5815286 ]))

Tout fonctionne bien.

Assemblage¶

On assemble une régression linéaire à coefficients positifs, une autre à coefficients négatifs, et une dernière pour assembler le tout. On commence par estimer une régression linéaire qu’on comparera au résultat final.

[16]:

df = load_wines_dataset()
df["color"] = df["color"].apply(lambda s: 1 if s == "red" else 0)
df.head()

[16]:

	fixed_acidity	volatile_acidity	citric_acid	residual_sugar	chlorides	free_sulfur_dioxide	total_sulfur_dioxide	density	pH	sulphates	alcohol	quality	color
0	7.4	0.70	0.00	1.9	0.076	11.0	34.0	0.9978	3.51	0.56	9.4	5	1
1	7.8	0.88	0.00	2.6	0.098	25.0	67.0	0.9968	3.20	0.68	9.8	5	1
2	7.8	0.76	0.04	2.3	0.092	15.0	54.0	0.9970	3.26	0.65	9.8	5	1
3	11.2	0.28	0.56	1.9	0.075	17.0	60.0	0.9980	3.16	0.58	9.8	6	1
4	7.4	0.70	0.00	1.9	0.076	11.0	34.0	0.9978	3.51	0.56	9.4	5	1

[17]:

X, y = df.drop("quality", axis=1), df["quality"]
X_train, X_test, y_train, y_test = train_test_split(X, y)

[18]:

lr = LinearRegression(fit_intercept=False)
lr.fit(X_train, y_train)
r2_score(y_test, lr.predict(X_test))

[18]:

0.3078620360903338

Si le modèle assemblé est équivalent, on devrait retrouver la même performance.

[19]:

stacking_lr = StackingRegressor(
    estimators=[
        ("plr", PositiveOrNegativeLinearRegression()),
        ("nlr", PositiveOrNegativeLinearRegression(positive=False)),
    ],
    final_estimator=LinearRegression(fit_intercept=False),
    cv=5,  # Nombre de folds pour la validation croisée
)
stacking_lr.fit(X_train, y_train)

[19]:

StackingRegressor(cv=5,
                  estimators=[('plr', PositiveOrNegativeLinearRegression()),
                              ('nlr',
                               PositiveOrNegativeLinearRegression(positive=False))],
                  final_estimator=LinearRegression(fit_intercept=False))

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

[20]:

r2_score(y_test, stacking_lr.predict(X_test))

[20]:

-5.418365357962336

Ca ne marche visiblement pas. Regardons les coefficients.

[21]:

for est in stacking_lr.estimators_:
    print(est, est.theta_)

PositiveOrNegativeLinearRegression() fixed_acidity            2.012826
volatile_acidity         0.093357
citric_acid              0.089653
residual_sugar           1.514515
chlorides                0.015391
free_sulfur_dioxide      8.550869
total_sulfur_dioxide    32.164952
density                  0.277868
pH                       0.899687
sulphates                0.148739
alcohol                  2.953378
color                    0.067369
dtype: float64
PositiveOrNegativeLinearRegression(positive=False) fixed_acidity           0.000000
volatile_acidity        0.000000
citric_acid             0.000000
residual_sugar          0.000000
chlorides              -0.016157
free_sulfur_dioxide     0.000000
total_sulfur_dioxide    0.000000
density                 0.000000
pH                      0.000000
sulphates               0.000000
alcohol                 0.000000
color                   0.000000
dtype: float64

L’algorithme d’apprentissage est instable, les coefficients sont parfois tous nuls. Il faudrait augmenter le nombre d’itérations.

Avec un pipeline¶

L’inconvénient de la méthode est que l’apprentissage est légèrement différent puisqu’il s’effectue avec une validation croisée. On peut supprimer ceci en créant son propre pipeline. Un pipeline empile des estimateurs avec pour seule constrainte que seul le dernier peut être un prédicteur, les autres sont des transformeurs.

Pour utiliser les prédictions d’un régresseur comme features d’un autre régression, il faut cacher le régresseur dans un transformeur.

[32]:

class PositiveOrNegativeLinearRegressionTransformer(TransformerMixin, BaseEstimator):
    def __init__(
        self, epsilon: float = 1.0, max_iter: int = 200, positive: bool = True
    ):
        super().__init__()
        self.epsilon = epsilon
        self.max_iter = max_iter
        self.positive = positive

    def fit(self, X, y):
        estimator = PositiveOrNegativeLinearRegression(
            epsilon=self.epsilon, max_iter=self.max_iter, positive=self.positive
        )
        estimator.fit(X, y)
        self.estimator_ = estimator
        self.n_targets_ = 1 if len(y.shape) == 1 else y.shape[0]
        return self

    def transform(self, X):
        return self.estimator_.predict(X).reshape((-1, self.n_targets_))


X = np.random.randn(5, 2)
y = X[:, 0] * 0.3 + X[:, 1] * 0.6 + np.random.randn(5) / 100
plr = PositiveOrNegativeLinearRegressionTransformer(max_iter=50, positive=True)
plr.fit(X, y)
plr.estimator_.theta_

[32]:

array([0.29527481, 0.60055273])

Tout fonctionne bien. Et maintenant un pipeline. Pour assembler les sorties de deux transformers, il faut utiliser un ColumnTransformer.

[51]:

model = Pipeline(
    [
        (
            "regposneg",
            ColumnTransformer(
                [
                    (
                        "pos",
                        PositiveOrNegativeLinearRegressionTransformer(),
                        make_column_selector(dtype_include=np.number),
                    ),
                    (
                        "neg",
                        PositiveOrNegativeLinearRegressionTransformer(positive=False),
                        make_column_selector(dtype_include=np.number),
                    ),
                ]
            ),
        ),
        ("norm", StandardScaler(with_mean=False)),
        ("lr", LinearRegression(fit_intercept=False)),
    ]
)

model.fit(X_train, y_train)
r2_score(y_test, model.predict(X_test))

[51]:

-4.884340735009395

C’est un peu mieux et surtout comparable. On regarde les coefficients, d’abord la régression à coefficients positifs.

[52]:

model.steps[0][1].transformers_[0][1].estimator_.theta_

[52]:

fixed_acidity            2.003196
volatile_acidity         0.092911
citric_acid              0.089224
residual_sugar           1.507269
chlorides                0.015317
free_sulfur_dioxide      8.509955
total_sulfur_dioxide    32.011052
density                  0.276539
pH                       0.895382
sulphates                0.148028
alcohol                  2.939247
color                    0.067047
dtype: float64

La régression à coefficients négatifs.

[53]:

model.steps[0][1].transformers_[1][1].estimator_.theta_

[53]:

fixed_acidity           0.000000
volatile_acidity        0.000000
citric_acid             0.000000
residual_sugar          0.000000
chlorides              -0.027225
free_sulfur_dioxide     0.000000
total_sulfur_dioxide    0.000000
density                 0.000000
pH                      0.000000
sulphates               0.000000
alcohol                 0.000000
color                   0.000000
dtype: float64

Puis la régression finale.

[54]:

model.steps[2][1].coef_

[54]:

array([ 1.54518281, -1.22963995])

Le score $R^2$ est très faible, difficile d’en tirer une conclusion. Il est aussi très inférieur au $R^2$ obtenu avec une régression linéaire simple. Les deux modèles ne sont donc pas équivalents.

Notebook on github