PRML Chapitre 4 Implémentation Python de la régression logistique bayésienne

Dans l'ajustement de courbe, nous le traitons comme Bayes et calculons la distribution de prédiction postérieure, mais on a l'impression qu'une telle chose n'est pas tellement faite en classification, donc cette fois nous allons faire une régression logistique, qui est souvent utilisée dans la classification, plus bayésienne. Je voudrais implémenter le code pour calculer la distribution de prédiction postérieure.

Retour logistique Bayes

Comme je l'ai écrit ci-dessus, je n'ai pas vu beaucoup ou pas de code qui calcule la distribution prédite dans la classification (peut-être juste quelques codes que j'ai vus). Afin de calculer la distribution prédite de manière bayésienne, la distribution postérieure des paramètres de poids doit être utilisée pour intégrer et éliminer les paramètres. Cependant, comme la régression logistique utilise la fonction sigmoïde logistique, il n'est pas possible d'intégrer analytiquement les paramètres. Par conséquent, la distribution prédite est approximativement obtenue en utilisant l'approximation de Laplace.

Retour logistique

Tout d'abord, considérons la régression logistique, qui classe deux classes. Soit $ \ phi $ le vecteur caractéristique d'un point $ x $ et que ce point appartienne à la classe $ C_1 $

p(C_1|\phi) = y(\phi) = \sigma({\bf w}^\intercal\phi)

Il est exprimé comme. cependant,\sigma(\cdot)Est une fonction sigmoïde logistique\sigma(x)={1\over1+\exp(-x)}Et. Puisque nous considérons ici la classification à deux classes, les classesC_1Si vous connaissez la probabilité d'appartenir à l'autre classeC_2La probabilité d'appartenir àp(C_2|\phi(x)) = 1 - p(C_1|\phi(x))Donné dans. Par conséquent, après cela, la classeC_1Seule la probabilité d'appartenance à est prise en compte.

Ensemble de données d'entraînement $ \ {\ phi_n, t_n \} _ {n = 1} ^ N $ as $ \ phi_n = \ phi (x_n) $, $ t_n \ in \ {0,1 \} $ Étant donné que la fonction de vraisemblance est modélisée sur la distribution de Bernoulli

\begin{align}
p({\bf t}|{\bf\Phi},{\bf w}) &= \prod_{n=1}^N{\rm Bern}(t_n|y_n)\\
&= \prod_{n=1}^N y_n^{t_n}(1 - y_n)^{1 - t_n}
\end{align}

Sera. Cependant, $ {\ bf \ Phi} $ est une matrice de planification dans laquelle les vecteurs horizontaux $ \ phi ^ {\ rm T} $ sont disposés verticalement. Comme toujours, définir la fonction de coût comme une fonction de vraisemblance logarithmique négative prend la forme d'une entropie croisée.

E({\bf w}) = -\ln p({\bf t}|{\bf\Phi},{\bf w}) = -\sum_{n=1}^N\{t_n\ln y_n + (1 - t_n)\ln(1 - y_n)\}

Les résultats de la différenciation du premier et du second ordre de cette fonction de coût pour le paramètre $ {\ bf w} $ sont

\begin{align}
\nabla E({\bf w}) &= \sum_{n=1}^N (y_n - t_n) {\bf\phi}_n = {\bf\Phi}^\intercal({\bf y} - {\bf t})\\
{\bf H} = \nabla\nabla E({\bf w}) &= \sum_{n=1}^N y_n(1 - y_n)\phi_n\phi_n^\intercal = {\bf\Phi}^\intercal{\bf R}{\bf\Phi}
\end{align}

Cependant, la matrice $ {\ bf R} $ est une matrice diagonale dont les éléments sont $ R_ {nn} = y_n (1-y_n) $. En utilisant ces résultats,

{\bf w}^{new} = {\bf w}^{old} - {\bf H}^{-1}\nabla E({\bf w})

Au fur et à mesure, la mise à jour est répétée jusqu'à ce que le paramètre $ {\ bf w} $ converge.

Approximation de Laplace

Distribution antérieure pour le paramètre $ {\ bf w} $ $ p ({\ bf w}) = \ mathcal {N} ({\ bf w} | {\ bf 0}, \ alpha ^ {-1} {\ bf I }) Introduisez $ et utilisez l'approximation de Laplace pour trouver une distribution postérieure approximative sous la forme d'une fonction gaussienne. Du théorème de Bayes

p({\bf w}|{\bf t},{\bf\Phi}) \propto p({\bf t}|{\bf\Phi},{\bf w})p({\bf w})

Donc, si vous prenez ce logarithme

\ln p({\bf w}|{\bf t},{\bf\Phi}) = -{\alpha\over2}{\bf w}^\intercal{\bf w} + \sum_{n=1}^N\{t_n\ln y_n + (1 - t_n)\ln(1 - y_n)\}

Sera. Pour calculer l'approximation gaussienne à partir de cela, maximisez d'abord l'équation ci-dessus pour $ {\ bf w} $, puis le paramètre $ {\ bf w} _ {MAP} $ et le différentiel du second ordre $ S_N ^ {-1 } En utilisant $, la distribution postérieure approximative

q({\bf w}) = \mathcal{N}({\bf w}|{\bf w}_{MAP}, S_N)

Est requis comme.

Distribution de prédiction ex post

La distribution de prédiction postérieure pour le nouveau point $ \ phi $ peut être obtenue en marginalisant la distribution a posteriori des paramètres.

\begin{align}
p(C_1|\phi,{\bf t},{\bf\Phi}) &= \int p(C_1|\phi,{\bf w})p({\bf w}|{\bf t},{\bf\Phi}){\rm d}{\bf w}\\
&\approx \int \sigma({\bf w}^\intercal\phi)\mathcal{N}({\bf w}|{\bf w}_{MAP},S_N){\rm d}{\bf w}\\
&\approx \sigma\left({\mu\over\sqrt{1 + \pi\sigma^2/8}}\right)
\end{align}

cependant,

\begin{align}
\mu &= {\bf w}_{MAP}^\intercal\phi\\
\sigma^2 &= \phi^\intercal S_N\phi
\end{align}

Et. Deux approximations sont utilisées: l'approximation de Laplace de la distribution postérieure du paramètre $ {\ bf w} $, et l'approximation de la fonction Probit de la fonction sigmoïde logistique. Voir PRML pour la transformation de formule détaillée.

la mise en oeuvre

paquet

En plus des habituels matplotlib et numpy, il utilise une bibliothèque Python standard appelée itertools.

import itertools
import matplotlib.pyplot as plt
import numpy as np

Matrice de conception: fonctions polygonales

Conversion d'un point bidimensionnel $ (x_1, x_2) $ en un vecteur de caractéristiques polymorphes quadratiques

\phi(x_1,x_2) =
\begin{bmatrix}
1\\
x_1\\
x_2\\
x_1^2\\
x_1x_2\\
x_2^2
\end{bmatrix}

Ce sera. Vous trouverez ci-dessous une classe qui convertit un certain nombre de points bidimensionnels en une matrice de conception. Par exemple, la matrice de planification des caractéristiques polymorphes quadratiques des points bidimensionnels $ (2,5) $ et $ (-3,1) $

{\bf\Phi}=
\begin{bmatrix}
1 & 2 & 5 & 2^2 & 2 \times 5 & 5^2\\
1 & 3 & -1 & 3^2 & 3 \times (-1) & (-1)^2
\end{bmatrix}

On dirait.

class PolynomialFeatures(object):

    def __init__(self, degree=2):
        self.degree = degree

    def transform(self, x):
        x_t = x.transpose()
        features = [np.ones(len(x))]
        for degree in xrange(1, self.degree + 1):
            for items in itertools.combinations_with_replacement(x_t, degree):
                features.append(reduce(lambda x, y: x * y, items))
        return np.array(features).transpose()

Retour logistique

Il n'est pas nécessaire de séparer les classes qui effectuent la régression logistique et la régression logistique de Bayes, mais cette fois nous les avons séparées pour faire la différence entre les deux.

class LogisticRegression(object):

    def __init__(self, iter_max, alpha=0):
        self.iter_max = iter_max
        self.alpha = alpha

    def _sigmoid(self, a):
        return np.divide(1, 1 + np.exp(-a))

    def fit(self, X, t):
        self.w = np.zeros(np.size(X, 1))
        for i in xrange(self.iter_max):
            w = np.copy(self.w)
            y = self.predict_proba(X)
            grad = X.T.dot(y - t) + self.alpha * w
            hessian = (X.T.dot(np.diag(y * (1 - y))).dot(X)
                       + self.alpha * np.identity(len(w)))
            try:
                self.w -= np.linalg.solve(hessian, grad)
            except np.linalg.LinAlgError:
                break
            if np.allclose(w, self.w):
                break
            if i == self.iter_max - 1:
                print "weight parameter w may not have converged"

    def predict(self, X):
        return (self.predict_proba(X) > 0.5).astype(np.int)

    def predict_proba(self, X):
        return self._sigmoid(X.dot(self.w))

Retour logistique Bayes

Lors de l'estimation de la probabilité postérieure maximale d'un paramètre, on obtient non seulement la valeur la plus fréquente, mais également une variance approximative. La régression logistique ci-dessus a également recherché la matrice de Hesse (la matrice de précision des paramètres approximatifs), mais je n'ai stocké cette matrice dans aucune variable car elle n'utilise que cette matrice lors de l'estimation. Cependant, la régression logistique bayésienne utilise également la variance des paramètres pour calculer la distribution prédite et la sauvegarde. Ensuite, la méthode predict_dist calcule la distribution prédite à l'aide de la matrice de distribution de ce paramètre.

class BayesianLogisticRegression(LogisticRegression):

    def __init__(self, iter_max, alpha=0.1):
        super(BayesianLogisticRegression, self).__init__(iter_max, alpha)

    def fit(self, X, t):
        super(BayesianLogisticRegression, self).fit(X, t)
        y = self.predict_proba(X)
        hessian = (X.T.dot(np.diag(y * (1 - y))).dot(X)
                   + self.alpha * np.identity(len(self.w)))
        self.w_var = np.linalg.inv(hessian)

    def predict_dist(self, X):
        mu_a = X.dot(self.w)
        var_a = np.sum(X.dot(self.w_var) * X, axis=1)
        return self._sigmoid(mu_a / np.sqrt(1 + np.pi * var_a / 8))

Code entier

import itertools
import matplotlib.pyplot as plt
import numpy as np


class PolynomialFeatures(object):

    def __init__(self, degree=2):
        self.degree = degree

    def transform(self, x):
        x_t = x.transpose()
        features = [np.ones(len(x))]
        for degree in xrange(1, self.degree + 1):
            for items in itertools.combinations_with_replacement(x_t, degree):
                features.append(reduce(lambda x, y: x * y, items))
        return np.array(features).transpose()


class LogisticRegression(object):

    def __init__(self, iter_max, alpha=0):
        self.iter_max = iter_max
        self.alpha = alpha

    def _sigmoid(self, a):
        return np.divide(1, 1 + np.exp(-a))

    def fit(self, X, t):
        self.w = np.zeros(np.size(X, 1))
        for i in xrange(self.iter_max):
            w = np.copy(self.w)
            y = self.predict_proba(X)
            grad = X.T.dot(y - t) + self.alpha * w
            hessian = (X.T.dot(np.diag(y * (1 - y))).dot(X)
                       + self.alpha * np.identity(len(w)))
            try:
                self.w -= np.linalg.inv(hessian).dot(grad)
            except np.linalg.LinAlgError:
                break
            if np.allclose(w, self.w):
                break
            if i == self.iter_max - 1:
                print "weight parameter w may not have converged"

    def predict(self, X):
        return (self.predict_proba(X) > 0.5).astype(np.int)

    def predict_proba(self, X):
        return self._sigmoid(X.dot(self.w))


class BayesianLogisticRegression(LogisticRegression):

    def __init__(self, iter_max, alpha=0.1):
        super(BayesianLogisticRegression, self).__init__(iter_max, alpha)

    def fit(self, X, t):
        super(BayesianLogisticRegression, self).fit(X, t)
        y = self.predict_proba(X)
        hessian = (X.T.dot(np.diag(y * (1 - y))).dot(X)
                   + self.alpha * np.identity(len(self.w)))
        self.w_var = np.linalg.inv(hessian)

    def predict_dist(self, X):
        mu_a = X.dot(self.w)
        var_a = np.sum(X.dot(self.w_var) * X, axis=1)
        return self._sigmoid(mu_a / np.sqrt(1 + np.pi * var_a / 8))


def create_data_set():
    x = np.random.normal(size=50).reshape(-1, 2)
    y = np.random.normal(size=50).reshape(-1, 2)
    y += np.array([2., 2.])
    return (np.concatenate([x, y]), np.concatenate([np.zeros(25), np.ones(25)]))


def main():
    x, labels = create_data_set()
    colors = ['blue', 'red']
    plt.scatter(x[:, 0], x[:, 1], c=[colors[int(label)] for label in labels])

    features = PolynomialFeatures(degree=3)

    classifier = BayesianLogisticRegression(iter_max=100, alpha=0.1)
    classifier.fit(features.transform(x), labels)

    X_test, Y_test = np.meshgrid(np.linspace(-2, 4, 100), np.linspace(-2, 4, 100))
    x_test = np.array([X_test, Y_test]).transpose(1, 2, 0).reshape(-1, 2)
    probs = classifier.predict_proba(features.transform(x_test))
    Probs = probs.reshape(100, 100)
    dists = classifier.predict_dist(features.transform(x_test))
    Dists = dists.reshape(100, 100)
    levels = np.linspace(0, 1, 5)
    cp = plt.contour(X_test, Y_test, Probs, levels, colors='k', linestyles='dashed')
    plt.clabel(cp, inline=True, fontsize=10)
    plt.contourf(X_test, Y_test, Dists, levels, alpha=0.5)
    plt.colorbar()
    plt.xlim(-2, 4)
    plt.ylim(-2, 4)
    plt.show()


if __name__ == '__main__':
    main()

résultat

• Données d'apprentissage: points bleus (classe 0), rouges (classe 1)
• ** Ligne pointillée noire: sortie de régression logistique **, une ligne avec une probabilité d'appartenir à la classe 1 à partir du haut 0,75,0,5,0,25, bien qu'il soit difficile à voir
• ** Codage couleur: sortie de régression logistique bayésienne **, chacun avec une probabilité d'appartenir à la classe 1
• Rouge: 0,75 ~ --Jaune: 0,5 à 0,75 --Mer: 0,25 ~ 0,5
• Bleu: ~ 0,25

En regardant cela, la ligne avec une probabilité de 0,5 (la ligne pointillée au milieu et la limite entre le bleu clair et le jaune) est commune dans les deux cas, tandis que les lignes de 0,25 et 0,75 sont plus de 0,5 pour ceux qui les traitent comme Bayes. Il est loin de la ligne et nous dit qu'on ne sait pas à quelle classe il appartient. En outre, lorsqu'il n'y a pas de points de données dans le coin supérieur gauche de la figure, la distribution prévue devient plus ambiguë que dans d'autres régions. D'autre part, comme il y a de nombreux points de données autour du centre, la largeur de 0,25 à 0,75 est étroite. Après avoir calculé la probabilité d'appartenir à la classe 1 par l'une ou l'autre méthode, elle est affectée à l'une ou l'autre des classes par un critère, mais si le critère est la minimisation des erreurs de classification, cela ne change pas celui qui est utilisé. Parce que la droite avec une probabilité de 0,5 est commune. Cependant, si vous souhaitez utiliser des critères plus complexes, vous voudrez peut-être le traiter comme bayésien.

À la fin

Cette fois, nous avons implémenté une méthode pour effectuer une régression logistique bayésienne en utilisant l'approximation de Laplace. L'avantage de le traiter de manière bayésienne est qu'il se comporte différemment dans les zones où il y a beaucoup de points de données et où il n'y en a pas. En plus d'utiliser l'approximation de Laplace, il semble y avoir une méthode pour dériver la distribution prévue en utilisant la variante de Bayes. Le chapitre 10 de PRML utilise également la variante de Bayes pour estimer le superparamètre $ \ alpha $. Je pense que je vais mettre cela en œuvre également.