Cette fois, j'ai implémenté l'approximation des preuves évoquée un peu à la fin de la première implémentation PRML. Cette méthode est également appelée Bayes empirique, le deuxième type d'estimation la plus probable, et ainsi de suite. Lorsque vous utilisez des algorithmes d'apprentissage automatique, vous devez souvent définir des super paramètres. En utilisant cette méthode, les super paramètres qui ont été décidés jusqu'à présent sont automatiquement décidés à partir des données.

Régression linéaire bayésienne

Avant d'entrer dans l'explication de l'approximation des preuves, parlons de la régression linéaire bayésienne.

Estimation la plus probable

Étant donné les données d'entraînement $ \ {x_i, t_i \} _ {i = 1} ^ N $, préparez le vecteur de caractéristiques $ {\ bf \ phi} $ (par exemple, $ {\ 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}\\

Comme, réduisez la fonction de vraisemblance suivante pour le paramètre $ {\ bf w} $. Cependant, $ {\ bf I} $ est la matrice unitaire de $ N \ fois N $, et $ \ beta $ est le paramètre de précision.

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

Distribution ex post

L'estimation la plus probable définit une distribution a priori dans le paramètre $ {\ bf w} $ car elle peut surcharger les données d'entraînement. Supposons donc que la distribution gaussienne de moyenne 0 et de variance $ \ alpha ^ {-1} $ soit une distribution antérieure, en utilisant le théorème de Bayes.

\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}

En tant que distribution postérieure pour le paramètre $ {\ bf w} $. cependant,

{\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}

Distribution de prédiction ex post

En utilisant la distribution a posteriori obtenue ci-dessus, la distribution prédite a posteriori de $ t $ pour la nouvelle entrée $ x $ est

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

Il peut être obtenu en intégrant et en éliminant le paramètre $ {\ bf w} $ comme dans. Cependant, la distribution est

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

En utilisant la distribution de prédiction postérieure obtenue de cette manière, la régression linéaire bayésienne était possible. Le comportement de cette régression dépend des super paramètres $ \ alpha, \ beta $. Si les valeurs des superparamètres que nous définissons sont trop grandes ou trop petites, les résultats de la régression peuvent ne pas être ceux que vous souhaitez.

Approximation des preuves

Dans ce qui précède, les super paramètres $ \ alpha, \ beta $ ont été définis comme des valeurs déterminées par l'homme, mais celles-ci sont estimées à partir des données. Lors de la recherche de la distribution postérieure du paramètre $ {\ 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)}

J'ai utilisé le théorème de Bayes comme. Lors de la recherche de la distribution postérieure, nous nous concentrons souvent uniquement sur la partie de la molécule du côté droit, mais en réalité, le dénominateur $ p ({\ bf t} | {\ bf \ Phi}, \ alpha, \ beta) $ est la fonction de preuve. Il est appelé, et il a la relation suivante.

(Fonction responsabilité)\times(Distribution antérieure)=(Fonction de preuve)\times(Distribution ex post)

Compte tenu de la signification de la fonction de preuve de $ p ({\ bf t} | {\ bf \ Phi}, \ alpha, \ beta) $, la fonction de preuve a les superparamètres donnés $ \ alpha, \ beta $. Montre à quel point il est facile de générer des données d'entraînement $ {\ bf t} $. Si la valeur de la fonction de preuve est petite, $ \ alpha, \ beta $ est difficile à générer les données, et inversement, si la valeur de la fonction de preuve est grande, les données sont faciles à générer. En d'autres termes, ** plus la valeur de la fonction de preuve est élevée, mieux c'est **. Approximativement en substituant les valeurs des super-paramètres qui maximisent la valeur de la fonction de preuve dans $ {\ hat \ alpha}, {\ hat \ beta} $ dans la formule de la distribution de prédiction postérieure.

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

Est l'approximation des preuves.

Maximiser les preuves

Trouvez $ \ alpha, \ beta $ qui maximise la fonction de preuve. Mettez à jour $ \ alpha, \ beta, \ gamma $ avec la valeur unique de $ {\ bf \ Phi} ^ {\ rm T} {\ bf \ Phi} $ comme $ \ lambda_i $ comme suit.

\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}

En répétant ces deux étapes, la fonction de preuve converge vers $ \ alpha = {\ hat \ alpha}, \ beta = {\ hat \ beta} $. Voir PRML pour la dérivation de cette formule.

la mise en oeuvre

Je réutilise le code que j'ai créé pour Bayes Curve Fitting.

Matrice de planification

La classe PolynomialFeatures transforme le vecteur d'entrée $ (x_1, \ cdots, x_N) ^ {\ rm T} $ en la matrice de planification $ {\ bf \ Phi} $ en définissant l'ordre.

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	Description de la méthode
__init__	Définir l'ordre des caractéristiques polymorphes
transform	Convertir l'entrée en matrice de conception

Régression linéaire bayésienne

La classe BayesianRegression estime la probabilité postérieure et calcule la distribution de prédiction postérieure comme décrit dans Bayesian Linear Regression.

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	Description de la méthode
__init__	\alpha,\betaDéfinissez la valeur de
fit	Paramètres{\bf w}Calculer la moyenne et la variance de la distribution postérieure de
predict	Calculer la valeur la plus fréquente de la distribution de prédiction postérieure
predict_dist	Calcul de la valeur la plus fréquente et de l'écart type de la distribution de prédiction postérieure

Approximation des preuves

La fonction de preuve est maximisée en mettant à jour alternativement le paramètre $ {\ bf w} $ et les super paramètres $ \ alpha, \ beta $. La mise à jour du paramètre $ {\ bf w} $ ne fait que la même chose que dans l'ajustement de la courbe bayésienne ci-dessus, nous réutilisons donc la méthode avec l'héritage de classe. La méthode de preuve calcule également la valeur de la fonction de preuve logarithmique. En comparant cette valeur, vous pouvez sélectionner automatiquement le meilleur modèle parmi plusieurs modèles.

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	Description de la méthode
__init__	\alpha,\betaDéfinissez les valeurs et le nombre maximum de fois pour les mettre à jour
fit	\alpha,\betaEstimation et paramètres{\bf w}Calculer la moyenne et la variance de la distribution postérieure de
evidence	Calculer la valeur de la fonction de preuve logarithmique

Code entier

Générez des points de données qui suivent le polynôme cubique $ x (x-5) (x + 5) $, puis effectuez un ajustement de courbe en utilisant l'approximation des preuves avec un modèle du polynôme d'ordre 0 au 7e pour obtenir les preuves du modèle. La distribution de prédiction postérieure est calculée à l'aide du modèle avec la plus grande valeur.

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

Définissez un polypole cubique $ x (x-5) (x + 5) $ (def func ())
Créez un point de données bruyant selon la fonction (x, t = create_toy_data (...))
Répétez les étapes 4 à 6 avec $ i = 0, \ cdots, 7 $
Créez une matrice de planification des entités polymorphes $ i $ order (X = features.transform (x))
Définissez les valeurs initiales de $ \ alpha et \ beta $ sur 100 (valeurs qui seraient inutiles si la fonction de preuve n'est pas maximisée), et effectuez un ajustement de courbe par approximation de preuve (regression.fit (X ,, X,) t) )
Calculez et enregistrez la valeur de la fonction de preuve logarithmique à l'étape 5 (ʻevidences.append (regression.evidence (X, t)) `)
Obtenez l'ordre du modèle avec la valeur la plus élevée de la fonction de preuve logarithmique (degree = np.argmax (evidences))
Obtenez le résultat de la régression dans cet ordre (regression = regressions [degree])
Les résultats ci-dessus sont illustrés

résultat

Comme le montre la figure 3.14 de PRML, la valeur de la fonction de preuve logarithmique lors de l'utilisation d'un polypole cubique comme modèle est la plus grande. De plus, la distribution prédite par le modèle du polynôme cubique est également intéressante. Puisque les valeurs initiales des super paramètres $ \ alpha et \ beta $ étaient toutes les deux fixées à 100, dans le cas d'un ajustement de courbe bayésienne ordinaire, $ \ alpha $ est trop grand et chaque composante de $ {\ bf w} $ La valeur est proche de 0, la variance de la distribution prédite est faible car $ \ beta $ est trop grande, et il aurait dû y avoir une distribution prédite postérieure qui ne pouvait pas du tout être ajustée aux données. Puisque nous avons estimé $ \ alpha, \ beta $ pour maximiser la fonction de preuve, la distribution de prédiction postérieure correspond également aux données.

À la fin

En utilisant cette approximation de preuves, il a estimé qu'il serait acceptable de définir un mauvais super paramètre comme valeur initiale. Avec la recherche de grille, je cherchais des super paramètres pour vérifier certains points de données d'entraînement, mais avec l'approximation des preuves, il est possible d'estimer des super paramètres avec un petit nombre d'essais en utilisant toutes les données d'entraînement. Je peux le faire. De plus, cette approximation des preuves est utilisée non seulement pour les modèles de régression linéaire, mais aussi pour estimer la régression logistique et les hyperparamètres des réseaux de neurones. Si j'ai une chance, je vais essayer de les mettre en œuvre.

PRML Chapitre 3 Preuve Implémentation approximative de Python

Régression linéaire bayésienne

Estimation la plus probable

Distribution ex post

Distribution de prédiction ex post

Approximation des preuves

Maximiser les preuves

la mise en oeuvre

Matrice de planification

Régression linéaire bayésienne

Approximation des preuves

Code entier

evidence_approximation.py

résultat

À la fin

`evidence_approximation.py`