Aperçu

Lorsque la distribution obtenue à partir de l'échantillon est multimodale, il n'est pas approprié de modéliser avec une distribution gaussienne simple. Les distributions multimodales peuvent être modélisées à l'aide d'une ** distribution gaussienne mixte ** qui combine plusieurs distributions gaussiennes. Cet article fournit un exemple d'utilisation de l'algorithme EM pour déterminer les paramètres d'une ** distribution gaussienne mixte **.

Premier à partir de la distribution de type pic unique

Estimation la plus probable

Soit l'échantillon $ x_n (n = 1,…, N) $. L'estimation la plus probable de la distribution gaussienne nous permet d'obtenir la moyenne et la variance sous la forme suivante.

\mu_{ML}=\frac{1}{N}\sum_{n=1}^N x_n \\\
\sigma^2_{ML}=\frac{1}{N-1}\sum_{n=1}^N (x_n-\mu_{ML})^2

La valeur de la variance utilise une estimation sans biais.

procédure

Générez 10000 échantillons à partir d'une distribution gaussienne avec une moyenne de 3,0 et une variance de 2,0.
Estimer la moyenne et la variance de l'échantillon obtenu en utilisant l'estimation la plus probable.

Code source

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
import math

#Fonction pour tracer l'histogramme
def draw_hist(xs, bins):
    plt.hist(xs, bins=bins, normed=True, alpha=0.5)

#Une fonction pour trouver la moyenne et la variance d'un échantillon donné en utilisant l'estimation la plus probable
def predict(data):
    mu = np.mean(data)
    var = np.var(data, ddof=1)  #Utilisez des estimations non biaisées
    return mu, var

def main():
    #Moyenne mu,Générer N échantillons qui suivent la distribution gaussienne de l'écart type std
    mu = 3.0
    v = 2.0
    std = math.sqrt(v)
    N = 10000
    data = np.random.normal(mu, std, N)
    #Effectuer l'estimation la plus probable,Trouvez la moyenne et la variance
    mu_predicted, var_predicted = predict(data)
    #Trouvez l'écart type à partir de la valeur de la variance
    std_predicted = math.sqrt(var_predicted)
    print("original: mu={0}, var={1}".format(mu, v))
    print(" predict: mu={0}, var={1}".format(mu_predicted, var_predicted))

    #Graphique des résultats
    draw_hist(data, bins=40)
    xs = np.linspace(min(data), max(data), 200)
    norm = mlab.normpdf(xs, mu_predicted, std_predicted)
    plt.plot(xs, norm, color="red")
    plt.xlim(min(xs), max(xs))
    plt.xlabel("x")
    plt.ylabel("Probability")
    plt.show()


if __name__ == '__main__':
    main()

Résultat d'exécution

L'histogramme (bleu) représente l'échantillon et la ligne rouge représente la distribution gaussienne modélisée à l'aide des valeurs estimées.

original: mu=3.0, var=2.0
 predict: mu=2.98719564872, var=2.00297779707

Nous avons pu trouver un modèle de distribution gaussienne appropriée.

Distribution bimodale

base de données

Old Faithful - Utilise des données de printemps intermittentes. Vous pouvez les télécharger à partir du lien ci-dessous. Old Faithful Le contenu de l'ensemble de données est le suivant.

1ère ligne: durée de la dernière éruption (minutes)
2ème rangée: temps d'attente jusqu'à la prochaine éruption (minutes)

Cette fois, nous n'utiliserons que la durée d'éruption la plus récente (première colonne) comme échantillon. À partir de cet échantillon, la distribution bimodale suivante a été obtenue.

Le sentiment de voir la distribution Il semble qu'une simple distribution gaussienne ne puisse pas être modélisée. Maintenant, modélisons en utilisant une distribution gaussienne mixte qui combine deux distributions gaussiennes. Lors de la modélisation, il est nécessaire de déterminer les paramètres de moyenne $ \ mu_1, \ mu_2 $, variance $ \ sigma_1 ^ 2, \ sigma_2 ^ 2 $, et probabilité de mélange $ \ pi $. En supposant que la distribution gaussienne est $ \ phi (x | \ mu, \ sigma ^ 2) $, la distribution gaussienne mixte est donnée par l'équation suivante.

y=(1-\pi)\phi(x|\mu_1, \sigma^2_1)+\pi\phi(x|\mu_2, \sigma^2_2)

Algorithme EM

La moyenne et la variance de la distribution gaussienne peuvent être obtenues en résolvant analytiquement la maximisation de la fonction de vraisemblance (PRML. Voir PRML /) §2.3.4). Cependant, il est difficile de résoudre analytiquement la maximisation de la fonction de vraisemblance de la distribution gaussienne mixte, de sorte que la maximisation est effectuée à l'aide de l'algorithme EM, qui est l'une des méthodes d'optimisation. L'algorithme EM est un algorithme qui répète deux étapes, l'étape E et l'étape M.

Étape E (Attente): Calculez le taux de charge.
Étape M (Maximisation): Calculer la moyenne pondérée et la variance.

Voir le §8.5 de Les éléments de l'apprentissage statistique pour le contenu détaillé de l'algorithme. La version anglaise du PDF peut être téléchargée gratuitement.

Code source

# -*- coding: utf-8 -*-

import numpy as np
import math
import random
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab

#Moyenne m,Distribution gaussienne de la variance v
def gaussian(x, m, v):
    p = math.exp(- pow(x - m, 2) / (2 * v)) / math.sqrt(2 * math.pi * v)
    return p

#Étape E
def e_step(xs, ms, vs, p):
    burden_rates = []
    for x in xs:
        d = (1 - p) * gaussian(x, ms[0], vs[0]) + p * gaussian(x, ms[1], vs[1])
        n = p * gaussian(x, ms[1], vs[1])
        burden_rate = n / d
        burden_rates.append(burden_rate)
    return burden_rates


#Étape M
def m_step(xs, burden_rates):
    d = sum([1 - r for r in burden_rates])
    n = sum([(1 - r) * x for x, r in zip(xs, burden_rates)])
    mu1 = n / d

    n = sum([(1 - r) * pow(x - mu1, 2) for x, r in zip(xs, burden_rates)])
    var1 = n / d

    d = sum(burden_rates)
    n = sum([r * x for x, r in zip(xs, burden_rates)])
    mu2 = n / d

    n = sum(r * pow(x - mu2, 2) for x, r in zip(xs, burden_rates))
    var2 = n / d

    N = len(xs)
    p = sum(burden_rates) / N

    return [mu1, mu2], [var1, var2], p


#Fonction de vraisemblance du journal
def calc_log_likelihood(xs, ms, vs, p):
    s = 0
    for x in xs:
        g1 = gaussian(x, ms[0], vs[0])
        g2 = gaussian(x, ms[1], vs[1])
        s += math.log((1 - p) * g1 + p * g2)
    return s

#Fonction pour tracer l'histogramme
def draw_hist(xs, bins):
    plt.hist(xs, bins=bins, normed=True, alpha=0.5)

def main():
    #Lire la première colonne de l'ensemble de données
    fp = open("faithful.txt")
    data = []
    for row in fp:
        data.append(float((row.split()[0])))
    fp.close()
    #mu, vs,Définir la valeur initiale de p
    p = 0.5
    ms = [random.choice(data), random.choice(data)]
    vs = [np.var(data), np.var(data)]
    T = 50  #Nombre d'itérations
    ls = []  #Enregistrer le résultat du calcul de la fonction de vraisemblance logarithmique
    #Algorithme EM
    for t in range(T):
        burden_rates = e_step(data, ms, vs, p)
        ms, vs, p = m_step(data, burden_rates)
        ls.append(calc_log_likelihood(data, ms, vs, p))

    print("predict: mu1={0}, mu2={1}, v1={2}, v2={3}, p={4}".format(
        ms[0], ms[1], vs[0], vs[1], p))
    #Graphique des résultats
    plt.subplot(211)
    xs = np.linspace(min(data), max(data), 200)
    norm1 = mlab.normpdf(xs, ms[0], math.sqrt(vs[0]))
    norm2 = mlab.normpdf(xs, ms[1], math.sqrt(vs[1]))
    draw_hist(data, 20)
    plt.plot(xs, (1 - p) * norm1 + p * norm2, color="red", lw=3)
    plt.xlim(min(data), max(data))
    plt.xlabel("x")
    plt.ylabel("Probability")

    plt.subplot(212)
    plt.plot(np.arange(len(ls)), ls)
    plt.xlabel("step")
    plt.ylabel("log_likelihood")
    plt.show()

if __name__ == '__main__':
    main()

Résultat d'exécution

predict: mu1=2.01860781706, mu2=4.27334342119, v1=0.0555176191851, v2=0.191024193785, p=0.651595365985

La figure ci-dessus montre les résultats de la modélisation d'un échantillon à partir d'un ensemble de données avec une distribution gaussienne mixte. La figure ci-dessous montre comment la probabilité logarithmique augmente lorsque l'algorithme EM est répété.

Autre

Valeur initiale dans l'algorithme EM

La moyenne était un échantillon choisi au hasard, la variance était la variance de l'échantillon et la probabilité de mélange était de 0,5 comme valeur initiale. Il existe plusieurs façons de sélectionner la valeur initiale, mais il semble que cela seul aboutira à un papier (?).

À propos de la distribution multi-pics

Puisque cet échantillon a une distribution bimodale, nous avons combiné deux distributions gaussiennes. Bien sûr, vous pouvez combiner trois distributions gaussiennes ou plus, mais le nombre de paramètres que vous devez déterminer augmentera.

À propos des dimensions de l'ensemble de données

Cette fois, nous avons utilisé un échantillon unidimensionnel, mais vous pouvez également utiliser l'algorithme EM pour les distributions gaussiennes multivariées. La distribution gaussienne multivariée a beaucoup plus de paramètres à déterminer que la distribution gaussienne unidimensionnelle (moyenne, matrice de covariance).

[PYTHON] Essayez de modéliser une distribution multimodale à l'aide de l'algorithme EM

Aperçu

Premier à partir de la distribution de type pic unique

Estimation la plus probable

procédure

Code source

Résultat d'exécution

Distribution bimodale

base de données

Algorithme EM

Code source

Résultat d'exécution

Autre

Valeur initiale dans l'algorithme EM

À propos de la distribution multi-pics

À propos des dimensions de l'ensemble de données