PRML Kapitel 2 Python-Implementierung von Student t-Distribution

Dieses Mal werden wir die wahrscheinlichste Schätzung der t-Verteilung des Schülers implementieren. Die t-Verteilung von Schülern ist als Verteilung bekannt, die für Ausreißer robuster ist als die Gauß-Verteilung. Wenn Sie sich jedoch gut erinnern, wurde diese Verteilung nie verwendet, daher ist sie eine gute Gelegenheit und ihre Robustheit. Ich werde es tatsächlich überprüfen. Zum Zeitpunkt der ersten PRML-Implementierung habe ich eine Regel festgelegt, dass außer der Standardbibliothek nur numpy möglich ist, diesmal jedoch eine Funktion, die nicht so vertraut ist wie die Digamma-Funktion. Ich habe auch scipy verwendet, weil es herauskam. Zum zweiten Mal in diesem Projekt werde ich so schnell wie möglich ein anderes Drittanbieterpaket als numpy verwenden und an die Zukunft erinnert.

Verteilung der Schüler

Die t-Verteilung des Schülers ist die Gaußsche Verteilung\mathcal{N}(x|\mu,\tau^{-1})Genauigkeitsparameter von\tauGammaverteilung als konjugierte Vorverteilung von{\rm Gam}(\tau|a,b)Es ist eine Verteilung, die durch Integrieren und Eliminieren der Genauigkeit unter Verwendung erhalten wird. Daraus kann interpretiert werden, dass die Student-t-Verteilung eine gemischte Gauß-Verteilung ist, bei der eine unendliche Anzahl von Gauß-Verteilungen mit demselben Mittelwert, aber unterschiedlichen Varianzen hinzugefügt werden. $p(x|\mu,a,b)=\int^{\infty}_0 \mathcal{N}(x|\mu,\tau^{-1}){\rm Gam}(\tau|a,b){\rm d}\tauDanach,\nu=2a,\lambda=a/bDann wird es in Form der t-Verteilung des Schülers sein. ${\rm St}(x|\mu,\lambda,\nu)={\Gamma(\nu/2+1/2)\over\Gamma(\nu/2)}\left({\lambda\over\pi\nu}\right)^{1/2}\left[1+{\lambda(x-\mu)^2\over\nu}\right]^{-\nu/2-1/2}$Nehmen Sie die wahrscheinlichste Schätzung vor, wenn einige Punkte x angegeben sind\mu,a,b(Oder\mu,\lambda,\nu)Ich würde gerne schätzen, aber im Gegensatz zu der wahrscheinlichsten Schätzung in der Gaußschen Verteilung ist sie nicht in geschlossener Form. Oberer, höher{\rm St}(x|\mu,\lambda,\nu)Es scheint, dass selbst wenn Sie die Parameter unterscheiden, es überhaupt keine schöne Form sein wird. In PRML wird der EM-Algorithmus verwendet, um die wahrscheinlichste t-Verteilung von Student abzuschätzen.(PRML-Übung 12.24)Es heißt zu verwenden. Der EM-Algorithmus ist eine häufig verwendete Methode zum Schätzen von Parametern in Situationen, in denen nicht beobachtete Daten vorliegen, und wird in Kapitel 9 von PRML vorgestellt. Präzisionsparameter\tau$Wenden Sie den EM-Algorithmus mit an. Die Berechnungsformel, die auf die t-Verteilung des Schülers angewendet wird, ist in PRML nicht enthalten, daher werde ich sie ein wenig berechnen. Ich würde mich freuen, wenn Sie auf Fehler hinweisen könnten.

E Schritt

Dies ist der E-Schritt (Expectation) des EM-Algorithmus. Legen Sie in diesem Schritt die Parameter ($ \ mu, a, b $) fest, die Sie schätzen möchten, und berechnen Sie aus welcher Verteilung (oder Genauigkeit $ \ tau $) die Gaußsche Verteilung, aus der jeder Abtastpunkt $ x_i $ abgetastet wird. Ich werde. $ p(\tau_i|x_i,\mu,a,b) = \mathcal{N}(x_i|\mu,\tau^{-1}){\rm Gam}(\tau_i|a,b)/const. $ Die Gammaverteilung wird für die vorherige Verteilung verwendet, und die Gauß-Funktion wird für die Wahrscheinlichkeitsfunktion verwendet. Wenn dies berechnet wird, wird die Gammaverteilung $ {\ rm Gam} (\ tau_i | a + {1 \ over2}, b + {1 \ over2} (x_i- \ mu) ^ 2) $ als hintere Verteilung erhalten. Aus dem erwarteten Wert dieser posterioren Verteilung kann der Genauigkeitsparameter für jeden Abtastpunkt als $ \ tau_i = {a + {1 \ over2} \ über b + {1 \ over2} (x_i- \ mu) ^ 2} $ erhalten werden. ..

M Schritt

Dies ist der M-Schritt (Maximierung) des EM-Algorithmus. Berechnen Sie die Protokollwahrscheinlichkeitsfunktion für die vollständigen Daten $ \ {x_i, \ tau_i \} $ und maximieren Sie sie für die Parameter $ \ mu, a, b $. Verwenden Sie für den Genauigkeitsparameter $ \ tau_i $ zu diesem Zeitpunkt den in Schritt E erhaltenen. $\sum_{i=1}^N\ln p(x_i,\tau_i|\mu,a,b) = \sum_{i=1}^N\ln\\{\mathcal{N}(x_i|\mu,\tau_i^{-1}){\rm Gam}(\tau_i|a,b)\\}Wenn Sie dies berechnen (ich bin nicht sicher) und den Teil extrahieren, an dem der zu schätzende Parameter beteiligt ist, $-{1\over2}\sum_i\tau_i(x_i-\mu)^2 + aN\ln b -N\ln\Gamma(a)+a\sum_i\ln\tau_i - b\sum_i\tau_i$ Unterscheiden Sie dies mit jedem Parameter und setzen Sie ihn auf 0, um die Gleichung zu lösen. $ \ mu = {\ sum_i \ tau_ix_i \ over \ sum_i \ tau_i} $ $ a = \ psi ^ {-1} (\ ln b + {1 \ over N} \ sum_i \ ln \ tau_i) $ Bei Differenzierung durch $ b = {aN \ over \ sum_i \ tau_i} $$ a kam die Digammafunktion $ \ psi (x) $ heraus. Wenn die Umkehrfunktion dieser Funktion tatsächlich existiert, könnte ich sie lösen, aber numpy und scipy haben nicht die Umkehrfunktion der Digammafunktion (die Digammafunktion ist in scipy). Daher wird Parameter a durch die Gradientenmethode nur wenig aktualisiert.

Wenn Sie sich zunächst die Lösungen von $ a und b $ ansehen, werden sie miteinander gemischt, sodass die obigen drei Gleichungen die logarithmische Wahrscheinlichkeitsfunktion wahrscheinlich nicht maximieren. Auf diese Weise wird der EM-Algorithmus, wenn die logarithmische Wahrscheinlichkeitsfunktion nicht vollständig maximiert ist, als verallgemeinerter EM-Algorithmus bezeichnet. Es scheint, dass der Grund, warum der Freiheitsgradparameter $ \ nu (= 2a) $ auf einen bestimmten Wert festgelegt ist, wenn die wahrscheinlichste Schätzung der t-Verteilung des Schülers durchgeführt wird, darin besteht, dass der M-Schritt genau ausgeführt wird. Wenn der Freiheitsgradparameter festgelegt ist, ist das verbleibende Schätzziel $ \ mu, \ lambda $, und ich denke, dass die logarithmische Wahrscheinlichkeitsfunktion leicht maximiert werden kann.

Ich werde einige Korrekturen und Ergänzungen basierend auf den Kommentaren unten vornehmen. Wenn es sich um den ursprünglichen E-Schritt handelt, den Genauigkeitsparameter\tau_iPosteriore Verteilung vonp(\tau_i|x_i,\mu,a,b) = {\rm Gam}(\tau_i|a+{1\over2},b+{1\over2}(x_i-\mu)^2)Endet mit dem berechneten Teil und im nachfolgenden M-Schritt mit der perfekten logarithmischen Wahrscheinlichkeitsfunktion\ln p(x_i,\tau_i|\mu,a,b)のそPosteriore Verteilung vonについての期待値\mathbb{E}[\sum_i\ln p(x_i,\tau_i|\mu,a,b)]Berechnen. Wenn das Berechnungsergebnis nur dort extrahiert wird, wo es sich auf den Parameter bezieht,-{1\over2}\sum_i\mathbb{E}\[\tau_i\](x_i-\mu)^2+a\sum_i\mathbb{E}\[\ln\tau_i\]+aN\lnb-b\sum_i\mathbb{E}[\tau_i]-N\ln\Gamma(a)UnddieFormderzumaximierendenFunktionistteilweiseunterschiedlich.DerEM-AlgorithmusimOriginalartikelisteineBeispielannäherungandieErwartungswertberechnung\mathbb{E}[\sum_i\lnp(x_i,\tau_i|\mu,a,b)]=\sum_i\lnp(x_i,\tau_i^{(sample)}|\mu,a,b)EsgibtjedochimmernureineStichprobengröße\tau_i^{(sample)}=\mathbb{E}[\tau_i]Ichwürdeesbegrüßen,wennSiedasBeispielinterpretierenkönnten.InderfolgendenAbbildungscheintesbiszueinemgewissenGradzufunktionieren,sodassdieseBeispielnäherungdieGenauigkeitmöglicherweisenichtsostarkbeeinflusst(ichbinwirklichfroh,wenndiespassiert).

Fluss der wahrscheinlichsten Schätzung

Um das Obige zusammenzufassen,

1. Stellen Sie die Anfangswerte der Parameter $ \ mu, a, b $ ein
2. Berechnen Sie den Genauigkeitsparameter $ \ tau $ für alle Stichprobenpunkte in Schritt E.
3. Aktualisieren Sie die Parameter $ \ mu, a, b $, sodass der Wert der logarithmischen Wahrscheinlichkeitsfunktion für vollständige Daten im Schritt M erhöht wird
4. Beenden Sie das Programm, wenn die Parameter konvergiert haben. Andernfalls kehren Sie zu Schritt E zurück

Verwenden Sie diesen verallgemeinerten EM-Algorithmus, um die Parameter der Student-t-Verteilung zu ermitteln.

Implementierung

import Importieren Sie Gamma- und Digammafunktionen aus matplotlib und numpy, um die Ergebnisse zu veranschaulichen, und zusätzlich scipy.

import matplotlib.pyplot as plt
import numpy as np
from scipy.special import gamma, digamma

Gaußsche Verteilung

Höchstwahrscheinlich Schätzung durch Gaußsche Verteilung nur zum Vergleich mit der Student-t-Verteilung

class Gaussian(object):

    def fit(self, x):
        self.mean = np.mean(x)
        self.var = np.var(x)

    def predict_proba(self, x):
        return (np.exp(-0.5 * (x - self.mean) ** 2 / self.var)
                / np.sqrt(2 * np.pi * self.var))

Verteilung der Schüler

Dies ist der Code, mit dem die t-Verteilung des Schülers am wahrscheinlichsten geschätzt wird. Die wahrscheinlichste Schätzung wird durch die Anpassungsmethode durchgeführt. Wiederholen Sie die Schritte E und M und beenden Sie den Vorgang, wenn die Parameter nicht mehr aktualisiert werden.

class StudentsT(object):

    def __init__(self, mean=0, a=1, b=1, learning_rate=0.01):
        self.mean = mean
        self.a = a
        self.b = b
        self.learning_rate = learning_rate

    def fit(self, x):
        while True:
            params = [self.mean, self.a, self.b]
            self._expectation(x)
            self._maximization(x)
            if np.allclose(params, [self.mean, self.a, self.b]):
                break

    def _expectation(self, x):
        self.precisions = (self.a + 0.5) / (self.b + 0.5 * (x - self.mean) ** 2)

    def _maximization(self, x):
        self.mean = np.sum(self.precisions * x) / np.sum(self.precisions)
        a = self.a
        b = self.b
        self.a = a + self.learning_rate * (
            len(x) * np.log(b)
            + np.log(np.prod(self.precisions))
            - len(x) * digamma(a))
        self.b = a * len(x) / np.sum(self.precisions)

    def predict_proba(self, x):
        return ((1 + (x - self.mean) ** 2/(2 * self.b)) ** (-self.a - 0.5)
                * gamma(self.a + 0.5)
                / (gamma(self.a) * np.sqrt(2 * np.pi * self.b)))

Ganzer Code

students_t_mle.py

import matplotlib.pyplot as plt
import numpy as np
from scipy.special import gamma, digamma


class Gaussian(object):

    def fit(self, x):
        self.mean = np.mean(x)
        self.var = np.var(x)

    def predict_proba(self, x):
        return (np.exp(-0.5 * (x - self.mean) ** 2 / self.var)
                / np.sqrt(2 * np.pi * self.var))


class StudentsT(object):

    def __init__(self, mean=0, a=1, b=1, learning_rate=0.01):
        self.mean = mean
        self.a = a
        self.b = b
        self.learning_rate = learning_rate

    def fit(self, x):
        while True:
            params = [self.mean, self.a, self.b]
            self._expectation(x)
            self._maximization(x)
            if np.allclose(params, [self.mean, self.a, self.b]):
                break

    def _expectation(self, x):
        self.precisions = (self.a + 0.5) / (self.b + 0.5 * (x - self.mean) ** 2)

    def _maximization(self, x):
        self.mean = np.sum(self.precisions * x) / np.sum(self.precisions)
        a = self.a
        b = self.b
        self.a = a + self.learning_rate * (
            len(x) * np.log(b)
            + np.log(np.prod(self.precisions))
            - len(x) * digamma(a))
        self.b = a * len(x) / np.sum(self.precisions)

    def predict_proba(self, x):
        return ((1 + (x - self.mean) ** 2/(2 * self.b)) ** (-self.a - 0.5)
                * gamma(self.a + 0.5)
                / (gamma(self.a) * np.sqrt(2 * np.pi * self.b)))


def main():

    # create toy data including outliers and plot histogram
    x = np.random.normal(size=20)
    x = np.concatenate([x, np.random.normal(loc=20., size=3)])
    plt.hist(x, bins=50, normed=1., label="samples")

    # prepare model
    students_t = StudentsT()
    gaussian = Gaussian()

    # maximum likelihood estimate
    students_t.fit(x)
    gaussian.fit(x)

    # plot results
    x = np.linspace(-5, 25, 1000)
    plt.plot(x, students_t.predict_proba(x), label="student's t", linewidth=2)
    plt.plot(x, gaussian.predict_proba(x), label="gaussian", linewidth=2)
    plt.legend()
    plt.show()


if __name__ == '__main__':
    main()

Ergebnis

Wenn Sie den obigen Code ausführen, erhalten Sie das folgende Ergebnis. Wie in Abbildung 2.16 von PRML gezeigt, ist die Anpassung an die t-Verteilung des Schülers sicherlich robust, selbst wenn es einen Ausreißer gibt. Der Durchschnitt der Student-t-Verteilung liegt bei etwa 0, aber der Durchschnitt der Gaußschen Verteilung liegt bei etwa 2,5, gezogen von den Ausreißern.

Am Ende

Es wurde bestätigt, dass die wahrscheinlichste Schätzung mit der t-Verteilung des Schülers robuster ist als die Gaußsche Verteilung als Modell. Wenn ich eine Chance habe, werde ich versuchen, das Problem der Kurvenregression mithilfe der t-Verteilung von Student zu lösen.