Dieses Mal habe ich die am Ende von [PRML first implement] erwähnte Evidenznäherung implementiert (http://qiita.com/cutting_the_Gordian_knot/items/555802600638f41b40c5). Diese Methode wird auch als empirische Bayes bezeichnet, die zweite Art der wahrscheinlichsten Schätzung und so weiter. Wenn Sie Algorithmen für maschinelles Lernen verwenden, müssen Sie häufig Superparameter einstellen. Mit dieser Methode werden die Superparameter, die bisher irgendwie festgelegt wurden, automatisch aus den Daten ermittelt.

Bayesianische lineare Regression

Bevor wir uns mit der Erklärung der Evidenznäherung befassen, wollen wir uns mit der Bayes'schen linearen Regression befassen.

Höchstwahrscheinlich Schätzung

Bereiten Sie unter Berücksichtigung der Trainingsdaten $ \ {x_i, t_i \} _ {i = 1} ^ N $ den Merkmalsvektor $ {\ bf \ phi} $ vor (z. B. $ {\ bf \ phi). } (x) = (1, x, x ^ 2, x ^ 3) ^ {\ rm T} $)

{\bf t} = 
\begin{bmatrix}
t_1\\
t_2\\
\vdots\\
t_N
\end{bmatrix}\\
{\bf\Phi} =
\begin{bmatrix}
\phi_1(x_1) & \phi_2(x_1) & \cdots & \phi_M(x_1)\\\
\phi_1(x_2) & \phi_2(x_2) & \cdots & \phi_M(x_2)\\\
\vdots & \vdots & \ddots & \vdots\\\
\phi_1(x_N) & \phi_2(x_N) & \cdots & \phi_M(x_N)
\end{bmatrix}\\
{\bf w} =
\begin{bmatrix}
w_1\\
w_2\\
\vdots\\
w_M
\end{bmatrix}\\

Minimieren Sie die folgende Wahrscheinlichkeitsfunktion für den Parameter $ {\ bf w} $. $ {\ Bf I} $ ist jedoch die Einheitsmatrix von $ N \ mal N $, und $ \ beta $ ist der Präzisionsparameter.

E({\bf w}) = \mathcal{N}({\bf t}|{\bf\Phi}{\bf w}, \beta^{-1}{\bf I})

Ex-post-Verteilung

Die wahrscheinlichste Schätzung legt eine vorherige Verteilung im Parameter $ {\ bf w} $ fest, da dadurch die Trainingsdaten möglicherweise überfordert werden. Nehmen wir also an, dass die Gaußsche Verteilung mit dem Mittelwert 0 und der Varianz $ \ alpha ^ {-1} $ die vorherige Verteilung unter Verwendung des Bayes-Theorems ist.

\begin{align}
p({\bf w}|{\bf\Phi},{\bf t}, \alpha, \beta) &= {\mathcal{N}({\bf t}|{\bf\Phi}{\bf w}, \beta^{-1}{\bf I}) \mathcal{N}({\bf w}|{\bf 0}, \alpha^{-1}{\bf I})\over p({\bf t}|{\bf\Phi}, \alpha, \beta)}\\
&= \mathcal{N}({\bf w}|{\bf m}_N, {\bf S}_N)
\end{align}

Als hintere Verteilung für den Parameter $ {\ bf w} $. Jedoch,

{\bf m}_N = \beta {\bf S}_N{\bf\Phi}^{\rm T}{\bf t}\\
{\bf S}_N^{-1} = \alpha{\bf I} + \beta{\bf\Phi}^{\rm T}{\bf\Phi}

Ex-post-Vorhersageverteilung

Unter Verwendung der oben erhaltenen posterioren Verteilung beträgt die posteriore vorhergesagte Verteilung von $ t $ für die neue Eingabe $ x $

p(t|x,{\bf\Phi},{\bf t},\alpha,\beta) = \int \mathcal{N}(t|{\bf w}^{\rm T}{\bf\phi}(x), \beta^{-1}) \mathcal{N}({\bf w}|{\bf m}_N, {\bf S}_N) {\rm d}{\bf w}\\
= \mathcal{N}(t|{\bf m}_N^{\rm T}{\bf\phi}(x), \sigma^2_N(x))

Sie können durch Integrieren und Eliminieren des Parameters $ {\ bf w} $ wie in erhalten werden. Die Verteilung ist jedoch

\sigma^2_N(x) = {1\over\beta} + \phi(x)^{\rm T}{\bf S}_N\phi(x)

Unter Verwendung der auf diese Weise erhaltenen posterioren Vorhersageverteilung war eine Bayes'sche lineare Regression möglich. Das Verhalten dieser Regression hängt von den Superparametern $ \ alpha, \ beta $ ab. Wenn die von uns festgelegten Superparameterwerte zu groß oder zu klein sind, entsprechen die Ergebnisse der Regression möglicherweise nicht Ihren Wünschen.

Evidenznäherung

Oben wurden die Superparameter $ \ alpha, \ beta $ als einige vom Menschen bestimmte Werte festgelegt, die jedoch aus den Daten geschätzt werden. Beim Finden der posterioren Verteilung des Parameters $ {\ bf w} $

p({\bf w}|{\bf\Phi},{\bf t}, \alpha, \beta) = {\mathcal{N}({\bf t}|{\bf\Phi}{\bf w}, \beta^{-1}{\bf I}) \mathcal{N}({\bf w}|{\bf 0}, \alpha^{-1}{\bf I})\over p({\bf t}|{\bf\Phi}, \alpha, \beta)}

Ich habe den Satz von Bayes als verwendet. Wenn wir die posteriore Verteilung finden, konzentrieren wir uns oft nur auf den Teil des Moleküls auf der rechten Seite, aber in Wirklichkeit ist der Nenner $ p ({\ bf t} | {\ bf \ Phi}, \ alpha, \ beta) $ die Beweisfunktion. Es heißt und hat die folgende Beziehung.

(Haftungsfunktion)\times(Vorherige Verteilung)=(Beweisfunktion)\times(Ex-post-Verteilung)

Unter Berücksichtigung der Bedeutung der Evidenzfunktion aus $ p ({\ bf t} | {\ bf \ Phi}, \ alpha, \ beta) $ hat die Evidenzfunktion die angegebenen Superparameter $ \ alpha, \ beta $. Zeigt, wie einfach es ist, Trainingsdaten $ {\ bf t} $ zu generieren. Wenn der Wert der Evidenzfunktion klein ist, ist es schwierig, die Daten mit $ \ alpha, \ beta $ zu generieren, und umgekehrt, wenn der Wert der Evidenzfunktion groß ist, sind die Daten leicht zu generieren. Mit anderen Worten, ** je höher der Wert der Evidenzfunktion ist, desto besser **. Ungefähr durch Ersetzen der Werte der Superparameter, die den Wert der Evidenzfunktion maximieren, in $ {\ hat \ alpha}, {\ hat \ beta} $ in der Formel der posterioren Vorhersageverteilung.

p(t|x,{\bf\Phi},{\bf t}) \approx p(t|x,{\bf\Phi},{\bf t},{\hat\alpha},{\hat\beta})

Ist die Evidenznäherung.

Beweise maximieren

Finden Sie $ \ alpha, \ beta $, das die Beweisfunktion maximiert. Aktualisieren Sie $ \ alpha, \ beta, \ gamma $ wie folgt mit dem eindeutigen Wert $ {\ bf \ Phi} ^ {\ rm T} {\ bf \ Phi} $ als $ \ lambda_i $.

\gamma=\sum_i{\beta\lambda_i\over\alpha+\beta\lambda_i}
\alpha={\gamma\over{\bf m}\_N^{\rm T}{\bf m}\_N}、\beta={N-\gamma\over\sum\_n\left\\{t\_n - {\rm m}\_N^{\rm T}{\bf\phi}(x\_n)\right\\}^2}

Durch Wiederholen dieser beiden Schritte konvergiert die Evidenzfunktion zu $ \ alpha = {\ hat \ alpha}, \ beta = {\ hat \ beta} $. Siehe PRML für die Ableitung dieser Formel.

Implementierung

Ich verwende den Code, den ich für [Bayes Curve Fitting] erstellt habe (http://qiita.com/cutting_the_Gordian_knot/items/555802600638f41b40c5).

Planungsmatrix

Die PolynomialFeatures-Klasse transformiert den Eingabevektor $ (x_1, \ cdots, x_N) ^ {\ rm T} $ in die Planungsmatrix $ {\ bf \ Phi} $, indem die Reihenfolge festgelegt wird.

import matplotlib.pyplot as plt
import numpy as np


class PolynomialFeatures(object):

    def __init__(self, degree):
        assert type(degree) is int, "%s is not int" % type(degree)
        self.degree = degree

    def transform(self, x):
        features = [x ** i for i in xrange(self.degree + 1)]
        return np.array(features).transpose()

PolynomialFeatures	Methodenbeschreibung
__init__	Legen Sie die Reihenfolge der polymorphen Merkmale fest
transform	Konvertieren Sie die Eingabe in eine Entwurfsmatrix

Bayesianische lineare Regression

Die BayesianRegression-Klasse schätzt die posteriore Wahrscheinlichkeit und berechnet die posteriore Vorhersageverteilung, wie in Bayesian Linear Regression beschrieben.

class BayesianRegression(object):

    def __init__(self, alpha=1., beta=1.):
        self.alpha = alpha
        self.beta = beta

    def fit(self, X, t):
        self.w_var = np.linalg.inv(
            self.alpha * np.identity(np.size(X, 1))
            + self.beta * X.T.dot(X))
        self.w_mean = self.beta * self.w_var.dot(X.T.dot(t))

    def predict(self, X):
        return X.dot(self.w_mean)

    def predict_dist(self, X):
        y = X.dot(self.w_mean)
        y_var = 1 / self.beta + np.sum(X.dot(self.w_var) * X, axis=1)
        y_std = np.sqrt(y_var)
        return y, y_std

BayesianRegression	Methodenbeschreibung
__init__	\alpha,\betaStellen Sie den Wert von ein
fit	Parameter{\bf w}Berechnen Sie den Mittelwert und die Varianz der posterioren Verteilung von
predict	Berechnen Sie den häufigsten Wert der posterioren Vorhersageverteilung
predict_dist	Berechnung des häufigsten Wertes und der Standardabweichung der posterioren Vorhersageverteilung

Evidenznäherung

Die Evidenzfunktion wird maximiert, indem abwechselnd der Parameter $ {\ bf w} $ und die Superparameter $ \ alpha, \ beta $ aktualisiert werden. Das Aktualisieren des Parameters $ {\ bf w} $ bewirkt nur dasselbe wie in der obigen Bayes'schen Kurvenanpassung. Daher verwenden wir die Methode mit Klassenvererbung wieder. Die Evidenzmethode berechnet auch den Wert der logarithmischen Evidenzfunktion. Durch Vergleichen dieses Werts können Sie automatisch das beste Modell aus mehreren Modellen auswählen.

class EvidenceApproximation(BayesianRegression):

    def __init__(self, iter_max=100, alpha=1., beta=1.):
        self.iter_max = iter_max
        self.alpha = alpha
        self.beta = beta

    def fit(self, X, t):
        M = X.T.dot(X)
        eigenvalues = np.linalg.eigvalsh(M)
        for i in xrange(self.iter_max):
            params = [self.alpha, self.beta]
            super(EvidenceApproximation, self).fit(X, t)
            self.gamma = np.sum(self.beta * eigenvalues / (self.alpha + self.beta * eigenvalues))
            self.alpha = self.gamma / self.w_mean.dot(self.w_mean)
            self.beta = (len(t) - self.gamma) / np.sum((t - X.dot(self.w_mean)) ** 2)
            if np.allclose(params, [self.alpha, self.beta]):
                break
        super(EvidenceApproximation, self).fit(X, t)

    def evidence(self, X, t):
        M = X.T.dot(X)
        return (len(M) * np.log(self.alpha)
                + len(t) * np.log(self.beta)
                - self.beta * np.sum((t - X.dot(self.w_mean)) ** 2)
                - np.linalg.slogdet(self.alpha + self.beta * M)[1])

EvidenceApproximation	Methodenbeschreibung
__init__	\alpha,\betaLegen Sie die Werte und die maximale Häufigkeit fest, mit der sie aktualisiert werden sollen
fit	\alpha,\betaSchätzung und Parameter{\bf w}Berechnen Sie den Mittelwert und die Varianz der posterioren Verteilung von
evidence	Berechnen Sie den Wert der logarithmischen Evidenzfunktion

Ganzer Code

Generieren Sie Datenpunkte, die dem kubischen Polynom $ x (x-5) (x + 5) $ folgen, und führen Sie dann eine Kurvenanpassung mithilfe der Evidenznäherung mit einem Modell des Polynoms 0. bis 7. Ordnung durch, um den Modellbeweis zu erhalten. Die posteriore Vorhersageverteilung wird unter Verwendung des Modells mit dem größten Wert berechnet.

`evidence_approximation.py`


import matplotlib.pyplot as plt
import numpy as np


class PolynomialFeatures(object):

    def __init__(self, degree):
        assert type(degree) is int, "%s is not int" % type(degree)
        self.degree = degree

    def transform(self, x):
        features = [x ** i for i in xrange(self.degree + 1)]
        return np.array(features).transpose()


class BayesianRegression(object):

    def __init__(self, alpha=1., beta=1.):
        self.alpha = alpha
        self.beta = beta

    def fit(self, X, t):
        self.w_var = np.linalg.inv(
            self.alpha * np.identity(np.size(X, 1))
            + self.beta * X.T.dot(X))
        self.w_mean = self.beta * self.w_var.dot(X.T.dot(t))

    def predict(self, X):
        return X.dot(self.w_mean)

    def predict_dist(self, X):
        y = X.dot(self.w_mean)
        y_var = 1 / self.beta + np.sum(X.dot(self.w_var) * X, axis=1)
        y_std = np.sqrt(y_var)
        return y, y_std


class EvidenceApproximation(BayesianRegression):

    def __init__(self, iter_max=100, alpha=1., beta=1.):
        self.iter_max = iter_max
        self.alpha = alpha
        self.beta = beta

    def fit(self, X, t):
        M = X.T.dot(X)
        eigenvalues = np.linalg.eigvalsh(M)
        for i in xrange(self.iter_max):
            params = [self.alpha, self.beta]
            super(EvidenceApproximation, self).fit(X, t)
            self.gamma = np.sum(self.beta * eigenvalues / (self.alpha + self.beta * eigenvalues))
            self.alpha = self.gamma / self.w_mean.dot(self.w_mean)
            self.beta = (len(t) - self.gamma) / np.sum((t - X.dot(self.w_mean)) ** 2)
            if np.allclose(params, [self.alpha, self.beta]):
                break
        super(EvidenceApproximation, self).fit(X, t)

    def evidence(self, X, t):
        M = X.T.dot(X)
        return (len(M) * np.log(self.alpha)
                + len(t) * np.log(self.beta)
                - self.beta * np.sum((t - X.dot(self.w_mean)) ** 2)
                - np.linalg.slogdet(self.alpha + self.beta * M)[1])


def create_toy_data(func, low=0, high=1, size=10, sigma=1.):
    x = np.random.uniform(low, high, size)
    t = func(x) + np.random.normal(scale=sigma, size=size)
    return x, t


def main():

    def func(x):
        return x * (x - 5) * (x + 5)

    x, t = create_toy_data(func, low=-5, high=5, size=20, sigma=5.)
    plt.scatter(x, t, s=50, alpha=0.5, label="observation")

    evidences = []
    regressions = []
    for i in xrange(8):
        features = PolynomialFeatures(degree=i)
        X = features.transform(x)
        regression = EvidenceApproximation(alpha=100., beta=100.)
        regression.fit(X, t)
        evidences.append(regression.evidence(X, t))
        regressions.append(regression)
    degree = np.argmax(evidences)
    regression = regressions[degree]

    x_test = np.linspace(-5, 5, 100)
    X_test = PolynomialFeatures(degree=int(degree)).transform(x_test)
    y, y_std = regression.predict_dist(X_test)

    plt.plot(x_test, func(x_test), color="blue", label="x(x-5)(x+5)")
    plt.plot(x_test, y, color="red", label="predict_mean")
    plt.fill_between(x_test, y - y_std, y + y_std, color="pink", alpha=0.5, label="predict_std")
    plt.legend()
    plt.show()

    plt.plot(evidences)
    plt.title("Model evidence")
    plt.xlabel("M")
    plt.ylabel("evidence")
    plt.show()


if __name__ == '__main__':
    main()

Definieren Sie ein kubisches Polypoly $ x (x-5) (x + 5) $ (def func ())
Erstellen Sie einen verrauschten Datenpunkt gemäß der Funktion (x, t = create_toy_data (...))
Wiederholen Sie die Schritte 4 bis 6 mit $ i = 0, \ cdots, 7 $
Erstellen Sie eine Planungsmatrix mit polymorphen Merkmalen der Ordnung $ i $ (X = features.transform (x)).
Setzen Sie die Anfangswerte von $ \ alpha und \ beta $ auf 100 (Werte, die nutzlos wären, wenn die Evidenzfunktion nicht maximiert wäre), und führen Sie eine Kurvenanpassung durch Evidenznäherung durch (regression.fit (X ,,) t) )
Berechnen und speichern Sie den Wert der logarithmischen Evidenzfunktion in Schritt 5 (evidences.append (Regression.evidence (X, t))).
Ermitteln Sie die Reihenfolge des Modells mit dem höchsten Wert der logarithmischen Evidenzfunktion (Grad = np.argmax (Beweise)).
Holen Sie sich das Ergebnis der Regression in dieser Reihenfolge (Regression = Regressionen [Grad])
Die obigen Ergebnisse sind dargestellt

Ergebnis

Wie in Abbildung 3.14 von PRML gezeigt, ist der Wert der logarithmischen Evidenzfunktion bei Verwendung eines kubischen Polypolys als Modell am größten. Auch die vorhergesagte Verteilung des kubischen Polynoms durch das Modell ist gut. Da die Anfangswerte der Superparameter $ \ alpha und \ beta $ beide auf 100 gesetzt wurden, ist $ \ alpha $ im Fall einer normalen Bayes'schen Kurvenanpassung zu groß und jede Komponente von $ {\ bf w} $ Der Wert liegt nahe bei 0, die Varianz der vorhergesagten Verteilung ist gering, da $ \ beta $ zu groß ist, und es hätte eine hintere vorhergesagte Verteilung geben müssen, die überhaupt nicht an die Daten angepasst werden konnte. Da wir $ \ alpha, \ beta $ geschätzt haben, um die Evidenzfunktion zu maximieren, passt die posteriore Vorhersageverteilung auch zu den Daten.

Am Ende

Unter Verwendung dieser Evidenznäherung schätzte er, dass es in Ordnung wäre, einen schlechten Superparameter als Anfangswert festzulegen. Bei der Rastersuche suchte ich nach Superparametern, um einige Punkte der Trainingsdaten zu überprüfen. Mit der Evidenznäherung ist es jedoch möglich, Superparameter mit einer kleinen Anzahl von Versuchen unter Verwendung aller Trainingsdaten zu schätzen. Ich kann es schaffen Darüber hinaus wird diese Evidenznäherung nicht nur für lineare Regressionsmodelle verwendet, sondern auch zur Schätzung der logistischen Regression und der Hyperparameter neuronaler Netze. Wenn ich eine Chance habe, werde ich sie auch umsetzen.

PRML Kapitel 3 Evidence Ungefähre Python-Implementierung

Bayesianische lineare Regression

Höchstwahrscheinlich Schätzung

Ex-post-Verteilung

Ex-post-Vorhersageverteilung

Evidenznäherung

Beweise maximieren

Implementierung

Planungsmatrix

Bayesianische lineare Regression

Evidenznäherung

Ganzer Code

evidence_approximation.py

Ergebnis

Am Ende

`evidence_approximation.py`