Cette fois, nous allons implémenter l'estimation la plus probable de la distribution t de Student. La distribution t de Student est bien connue comme une distribution plus robuste aux valeurs aberrantes que la distribution gaussienne, mais si vous vous souvenez bien, cette distribution n'a jamais été utilisée, c'est donc une bonne opportunité et sa robustesse. Je vais vérifier. Au moment de la Première implémentation PRML, j'ai défini une règle selon laquelle seul numpy est possible sauf pour la bibliothèque standard, mais cette fois, il s'agit d'une fonction inconnue appelée fonction digamma. J'ai aussi utilisé scipy parce qu'il est sorti. Pour la deuxième fois de ce projet, j'utiliserai un package tiers autre que numpy dès que possible, et je me rappellerai de l'avenir.

Distribution t de Student

La distribution t de Student est la distribution gaussienne\mathcal{N}(x|\mu,\tau^{-1})Paramètre de précision de\tauDistribution gamma en tant que distribution a priori conjuguée de{\rm Gam}(\tau|a,b)C'est une distribution obtenue en intégrant et en éliminant la précision à l'aide. À partir de là, on peut interpréter que la distribution t de Student est une distribution gaussienne mixte dans laquelle un nombre infini de distributions gaussiennes avec la même moyenne mais des variances différentes sont ajoutés. $p(x|\mu,a,b)=\int^{\infty}_0 \mathcal{N}(x|\mu,\tau^{-1}){\rm Gam}(\tau|a,b){\rm d}\tauAprès ça,\nu=2a,\lambda=a/bEnsuite, ce sera sous la forme de la distribution t de l'étudiant. ${\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}$Faire l'estimation la plus probable lorsque certains points x sont donnés\mu,a,b(Ou\mu,\lambda,\nu)Je voudrais estimer, mais contrairement au cas de l'estimation la plus probable avec une distribution gaussienne, elle n'est pas sous forme fermée. plus haut{\rm St}(x|\mu,\lambda,\nu)Il semble que même si vous différenciez les paramètres, ce ne sera pas du tout une belle forme. Dans PRML, l'algorithme EM est utilisé pour estimer la distribution t la plus probable de Student.(Exercice PRML 12.24)Il dit à utiliser. L'algorithme EM est une technique couramment utilisée pour estimer les paramètres dans des situations où il y a des données non observées et est présenté au chapitre 9 du PRML. Paramètres de précision\tau$Appliquez l'algorithme EM avec. La formule de calcul appliquée à la distribution t de l'étudiant n'est pas incluse dans PRML, je vais donc la calculer un peu. Je vous serais reconnaissant si vous pouviez signaler des erreurs.

Étape E

Il s'agit de l'étape E (Attente) de l'algorithme EM. Dans cette étape, fixez les paramètres que vous souhaitez estimer ($ \ mu, a, b $) et calculez à partir de quelle variance (ou précision $ \ tau $) la distribution gaussienne à partir de laquelle chaque point d'échantillonnage $ x_i $ est échantillonné. Je vais. $ p(\tau_i|x_i,\mu,a,b) = \mathcal{N}(x_i|\mu,\tau^{-1}){\rm Gam}(\tau_i|a,b)/const. $ La distribution gamma est utilisée pour la distribution antérieure et la fonction de Gauss est utilisée pour la fonction de vraisemblance. Lorsque cela est calculé, la distribution gamma $ {\ rm Gam} (\ tau_i | a + {1 \ over2}, b + {1 \ over2} (x_i- \ mu) ^ 2) $ est obtenue comme distribution postérieure. À partir de la valeur attendue de cette distribution postérieure, le paramètre de précision pour chaque point d'échantillonnage peut être obtenu comme $ \ tau_i = {a + {1 \ over2} \ over b + {1 \ over2} (x_i- \ mu) ^ 2} $. ..

Étape M

Il s'agit de l'étape M (Maximisation) de l'algorithme EM. Calculez la fonction de vraisemblance du journal pour les données complètes $ \ {x_i, \ tau_i \} $ et maximisez-la pour les paramètres $ \ mu, a, b $. Pour le paramètre de précision $ \ tau_i $ à ce moment, utilisez celui obtenu à l'étape E. $\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)\\}Si vous calculez cela (je ne suis pas sûr) et extrayez la partie où le paramètre que vous souhaitez estimer est impliqué, $-{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$ Différenciez-le avec chaque paramètre et définissez-le comme égal à 0 pour résoudre l'équation. $ \ mu = {\ sum_i \ tau_ix_i \ over \ sum_i \ tau_i} $ $ a = \ psi ^ {-1} (\ ln b + {1 \ over N} \ sum_i \ ln \ tau_i) $ Lorsqu'elle est différenciée par $ b = {aN \ over \ sum_i \ tau_i} $ a, la fonction digamma $ psi (x) $ est obtenue. Si la fonction inverse de cette fonction existe réellement, je pourrais la résoudre, mais numpy et scipy n'ont pas la fonction inverse de la fonction digamma (la fonction digamma est en scipy). Par conséquent, le paramètre a sera peu mis à jour par la méthode du gradient.

Si vous regardez les solutions de $ a et b $ en premier lieu, elles sont mélangées l'une avec l'autre, donc les trois équations ci-dessus ne maximisent probablement pas la fonction de vraisemblance logarithmique. De cette manière, l'algorithme EM lorsque la fonction de vraisemblance logarithmique n'est pas complètement maximisée est appelé l'algorithme EM généralisé. Il semble que la raison pour laquelle le paramètre de degré de liberté $ \ nu (= 2a) $ est fixé à une certaine valeur lors de l'exécution de l'estimation la plus probable de la distribution t de l'étudiant est que l'étape M est exécutée exactement. Si le paramètre de degré de liberté est fixe, la cible d'estimation restante sera $ \ mu, \ lambda $, et je pense que la fonction de vraisemblance logarithmique peut être facilement maximisée.

Je ferai quelques corrections et suppléments en fonction des commentaires ci-dessous. S'il s'agit de l'étape E d'origine, le paramètre de précision\tau_iDistribution postérieure dep(\tau_i|x_i,\mu,a,b) = {\rm Gam}(\tau_i|a+{1\over2},b+{1\over2}(x_i-\mu)^2)Se termine par la partie calculée, et à l'étape M suivante, la fonction de vraisemblance logarithmique parfaite\ln p(x_i,\tau_i|\mu,a,b)のそDistribution postérieure deについての期待値\mathbb{E}[\sum_i\ln p(x_i,\tau_i|\mu,a,b)]Calculer. Si le résultat du calcul est extrait uniquement là où il est lié au paramètre,-{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),Etlaformedelafonctionàmaximiserestpartiellementdifférente.L'algorithmeEMdel'articleoriginalestunexempled'approximationducalculdelavaleurattendue\mathbb{E}[\sum_i\lnp(x_i,\tau_i|\mu,a,b)]=\sum_i\lnp(x_i,\tau_i^{(sample)}|\mu,a,b)Cependant,iln'yatoujoursqu'uneseuletailled'échantillonchacun\tau_i^{(sample)}=\mathbb{E}[\tau_i]J'apprécieraissivouspouviezinterprétercelaenutilisantl'échantillon.Danslafigureci-dessous,celasemblefonctionnerdansunecertainemesure,donccetexempled'approximationpeutnepasaffecterautantlaprécision(jesuisvraimentheureuxsicelaseproduit).

Flux de l'estimation la plus probable

Pour résumer ce qui précède,

Définissez les valeurs initiales des paramètres $ \ mu, a, b $
Calculez le paramètre de précision $ \ tau $ pour tous les points d'échantillonnage à l'étape E
Mettez à jour les paramètres $ \ mu, a, b $ afin que la valeur de la fonction de vraisemblance logarithmique pour des données complètes soit augmentée à l'étape M
Quittez si les paramètres ont convergé, sinon retournez à l'étape E

Utilisez cet algorithme EM généralisé pour trouver les paramètres de la distribution t de Student.

la mise en oeuvre

import Importez les fonctions gamma et digamma de matplotlib et numpy pour illustrer les résultats, et en plus scipy.

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

Distribution gaussienne

Estimation la plus probable par distribution gaussienne uniquement pour comparaison avec la distribution t de Student

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

Distribution t de Student

C'est le code pour estimer le plus probable par la distribution t de l'étudiant. L'estimation la plus probable est effectuée par la méthode d'ajustement. Répétez les étapes E et M et terminez lorsque les paramètres ne sont plus mis à jour.

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

Code entier

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

résultat

Si vous exécutez le code ci-dessus, vous obtiendrez le résultat suivant. Comme le montre la figure 2.16 de PRML, l'ajustement par la distribution t de Student est certainement robuste même s'il existe une valeur aberrante. La moyenne de la distribution t de Student est d'environ 0, mais la moyenne de la distribution gaussienne est d'environ 2,5, tirée par les valeurs aberrantes.

À la fin

Il a été confirmé que l'estimation la plus probable avec la distribution t de Student est plus robuste que la distribution gaussienne comme modèle. Si j'ai une chance, j'essaierai de résoudre le problème de régression de courbe en utilisant la distribution t de Student.

PRML Chapter 2 Student t-Distribution Python Implementation

Distribution t de Student

Étape E

Étape M

Flux de l'estimation la plus probable

la mise en oeuvre

Distribution gaussienne

Distribution t de Student

Code entier

students_t_mle.py

résultat

À la fin

`students_t_mle.py`