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.

students_t.png

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

result.png 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.

Recommended Posts

Régression linéaire avec distribution t de Student
Régression linéaire avec statsmodels
Régression avec un modèle linéaire
[Python] Régression linéaire avec scicit-learn
Régression linéaire robuste avec scikit-learn
Régression linéaire
PRML Chapter 2 Student t-Distribution Python Implementation
Prédire l'été chaud avec un modèle de régression linéaire
Régression linéaire d'apprentissage automatique
Introduction à l'hypothèse Tensorflow-About et au coût de la régression linéaire
Programmation linéaire avec PuLP
Effectuer une analyse de régression avec NumPy
Essayez la régression avec TensorFlow
Essayez d'implémenter la régression linéaire à l'aide de Pytorch avec Google Colaboratory
Régression du noyau avec Numpy uniquement
Machine Learning: Supervision - Régression linéaire
Analyse de régression multiple avec Keras
3. Distribution normale avec un réseau neuronal!
Ridge retour avec Mllib à Pyspark
Méthode de régression linéaire utilisant Numpy
Régression linéaire en ligne en Python
Implémentation de la régression logistique avec NumPy
Introduction à la modélisation statistique bayésienne avec python ~ Essai de régression linéaire avec MCMC ~