[PYTHON] Lineare Regression mit Student's t-Verteilung

Als ich wahrscheinlichste Schätzung des Schülers mit t-Verteilung implementierte, sagte ich, dass ich diese Verteilung auch auf die lineare Regression anwenden würde. Ich habe es diesmal versucht.

Lineare Regression

Unter der Annahme, dass $ \ {x_n, t_n \} _ {n = 1} ^ N $ die Trainingsdaten sind, wird bei linearer Regression $ {\ bf \ phi} (\ cdot) $ als Merkmalsvektor verwendet, und der quadratische Fehler ist Schätzen Sie den Mindestparameter $ {\ bf w} ^ * $.

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

Dies entspricht der Schätzung der Parameter, die die Wahrscheinlichkeitsfunktion maximieren, unter Verwendung der folgenden Gaußschen Verteilung, wenn sie unter Verwendung der Wahrscheinlichkeitsverteilung neu interpretiert werden.

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

Die Gaußsche Verteilung wird häufig auch in diesem Sinne verwendet, weist jedoch auch einige Nachteile auf, z. B. eine relativ unsichere Rauschunterdrückung. Dieses Mal schätzen wir die Parameter unter Verwendung der Student-t-Verteilung in der Wahrscheinlichkeitsfunktion.

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

Der Umriss der t-Verteilung des Schülers ist in der folgenden Abbildung dargestellt (aus PRML Abbildung 2.15). Wenn es $ \ nu \ to \ infty $ (grün) ist, hat es die gleiche Form wie die Gaußsche Verteilung. Im Vergleich zur grünen Verteilung haben die blauen und roten Verteilungen einen breiteren Saum, wodurch sie robuster gegen Ausreißer sind.

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()

Ergebnis

Im Modell der blauen Linie, das die Gaußsche Verteilung für die Wahrscheinlichkeitsfunktion verwendet, wird sie zu den Ausreißern gezogen, aber in der Schüler-t-Verteilung ist die lineare Regression robust, selbst wenn es Ausreißer gibt.