Régression linéaire avec distribution t de Student

Quand j'ai implémenté L'estimation la plus probable de Student avec la distribution t, j'ai dit que j'appliquerais également cette distribution à la régression linéaire. J'ai essayé cette fois.

Régression linéaire

En supposant que $ \ {x_n, t_n \} _ {n = 1} ^ N $ sont les données d'apprentissage, en régression linéaire, $ {\ bf \ phi} (\ cdot) $ est utilisé comme vecteur de caractéristiques, et l'erreur carrée est Estimez le paramètre minimum $ {\ bf w} ^ * $.

{\bf w}^* = \arg \min_{\bf w}\sum_{n=1}^N(t_n - {\bf w}^\top{\bf \phi}(x_n))^2

Cela revient à estimer les paramètres qui maximisent la fonction de vraisemblance en utilisant la distribution gaussienne suivante lorsqu'ils sont réinterprétés à l'aide de la distribution de probabilité.

\prod_{n=1}^N\mathcal{N}(t_n|{\bf w}^\top{\bf\phi}(x_n), \sigma^2)

La distribution gaussienne est également souvent utilisée dans ce sens, mais elle présente également certains inconvénients, tels que le fait d'être relativement peu sûr contre le bruit. Cette fois, nous estimons les paramètres en utilisant la distribution t de Student dans la fonction de vraisemblance.

\prod_{n=1}^N{\rm St}(t_n|{\bf w}^\top{\bf\phi}(x_n),\lambda,\nu)

Le contour de la distribution t de l'étudiant est montré dans la figure ci-dessous (à partir de la figure 2.15 de PRML). Quand il est $ \ nu \ to \ infty $ (vert), il a la même forme que la distribution gaussienne. Par rapport à la distribution verte, les distributions bleue et rouge ont un ourlet plus large, ce qui les rend plus robustes contre les valeurs aberrantes.

code

students_t_regression.py

import functools
import itertools
import matplotlib.pyplot as plt
import numpy as np
import scipy.special as sp


class PolynomialFeatures(object):

    def __init__(self, degree=2):
        assert isinstance(degree, int)
        self.degree = degree

    def transform(self, x):
        if x.ndim == 1:
            x = x[:, None]
        x_t = x.transpose()
        features = [np.ones(len(x))]
        for degree in range(1, self.degree + 1):
            for items in itertools.combinations_with_replacement(x_t, degree):
                features.append(functools.reduce(lambda x, y: x * y, items))
        return np.asarray(features).transpose()


class LinearRegressor(object):

    def fit(self, X, t):
        self.w = np.linalg.pinv(X) @ t

    def predict(self, X):
        return X @ self.w


class RobustRegressor(LinearRegressor):

    def __init__(self, precision=1., dof=1.):
        self.precision = precision
        self.dof = dof

    def fit(self, X, t, learning_rate=0.01):
        super().fit(X, t)
        params = np.hstack((self.w.ravel(), self.precision, self.dof))
        while True:
            E_precision, E_ln_precision = self._expect(X, t)
            self._maximize(X, t, E_precision, E_ln_precision, learning_rate)
            params_new = np.hstack((self.w.ravel(), self.precision, self.dof))
            if np.allclose(params, params_new):
                break
            params = np.copy(params_new)

    def _expect(self, X, t):
        sq_diff = (t - X @ self.w) ** 2
        E_eta = (self.dof + 1) / (self.dof + self.precision * sq_diff)
        E_ln_eta = (
            sp.digamma(0.5 * (self.dof + 1))
            - np.log(0.5 * (self.dof + self.precision * sq_diff)))
        return E_eta, E_ln_eta

    def _maximize(self, X, t, E_eta, E_ln_eta, learning_rate):
        sq_diff = (t - X @ self.w) ** 2
        self.w = np.linalg.inv(X.T @ (E_eta[:, None] * X)) @ X.T @ (E_eta * t)
        self.precision = 1 / np.mean(E_eta * sq_diff)
        N = len(t)
        self.dof += learning_rate * (
            N * np.log(0.5 * self.dof) + N
            - N * sp.digamma(0.5 * self.dof)
            + np.sum(E_ln_eta - E_eta))


def create_toy_data(n):
    x = np.linspace(-1, 1, n)
    y = np.random.randint(-5, 6) * x + np.random.normal(scale=0.1, size=n)
    y[np.random.randint(len(y))] += np.random.uniform(-10, 10)
    return x, y


def main():
    x_train, y_train = create_toy_data(10)

    feature = PolynomialFeatures(degree=1)
    X_train = feature.transform(x_train)

    linear_regressor = LinearRegressor()
    robust_regressor = RobustRegressor()
    linear_regressor.fit(X_train, y_train)
    robust_regressor.fit(X_train, y_train)

    x = np.linspace(-1, 1, 100)
    X = feature.transform(x)
    y_linear = linear_regressor.predict(X)
    y_robust = robust_regressor.predict(X)

    plt.scatter(
        x_train, y_train,
        facecolors="none", edgecolors="g", label="training data")
    plt.plot(x, y_linear, c="b", label="Gaussian")
    plt.plot(x, y_robust, c="r", label="Student's t")
    plt.legend(loc="best")
    plt.show()


if __name__ == '__main__':
    main()

résultat

Dans le modèle de ligne bleue qui utilise la distribution gaussienne pour la fonction de vraisemblance, elle est tirée vers les valeurs aberrantes, mais dans la distribution t de Student, même s'il existe des valeurs aberrantes, la régression linéaire est robuste.