Überblick

Wenn die von der Probe erhaltene Verteilung multimodal ist, ist es nicht angebracht, mit einer einfachen Gaußschen Verteilung zu modellieren. Multimodale Verteilungen können unter Verwendung einer ** gemischten Gaußschen Verteilung ** modelliert werden, die mehrere Gaußsche Verteilungen kombiniert. Dieser Artikel enthält ein Beispiel für die Verwendung des EM-Algorithmus zur Bestimmung der Parameter einer ** gemischten Gaußschen Verteilung **.

Zunächst aus der Einzelpeak-Typverteilung

Höchstwahrscheinlich Schätzung

Die Stichprobe sei $ x_n (n = 1,…, N) $. Die wahrscheinlichste Schätzung der Gaußschen Verteilung ermöglicht es uns, den Mittelwert und die Varianz in der folgenden Form zu erhalten.

\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

Der Wert der Varianz verwendet eine unvoreingenommene Schätzung.

Verfahren

Generieren Sie 10000 Proben aus einer Gaußschen Verteilung mit einem Mittelwert von 3,0 und einer Varianz von 2,0.
Schätzen Sie den Mittelwert und die Varianz aus der erhaltenen Stichprobe mit der wahrscheinlichsten Schätzung.

Quellcode

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

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

#Funktion zum Zeichnen des Histogramms
def draw_hist(xs, bins):
    plt.hist(xs, bins=bins, normed=True, alpha=0.5)

#Eine Funktion, die den Mittelwert und die Varianz einer bestimmten Stichprobe anhand der wahrscheinlichsten Schätzung ermittelt
def predict(data):
    mu = np.mean(data)
    var = np.var(data, ddof=1)  #Verwenden Sie unvoreingenommene Schätzungen
    return mu, var

def main():
    #Durchschnitt mu,Generieren Sie N Stichproben, die der Gaußschen Verteilung der Standardabweichung std folgen
    mu = 3.0
    v = 2.0
    std = math.sqrt(v)
    N = 10000
    data = np.random.normal(mu, std, N)
    #Führen Sie die wahrscheinlichste Schätzung durch,Finden Sie den Mittelwert und die Varianz
    mu_predicted, var_predicted = predict(data)
    #Finden Sie die Standardabweichung vom Wert der Varianz
    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))

    #Ergebnisplot
    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()

Ausführungsergebnis

Das Histogramm (blau) repräsentiert die Probe und die rote Linie repräsentiert die Gaußsche Verteilung, die unter Verwendung der geschätzten Werte modelliert wurde.

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

Wir konnten ein Modell einer geeigneten Gaußschen Verteilung finden.

Bimodale Verteilung

Datensatz

Old Faithful-Uses intermittierende Frühlingsdaten. Sie können sie über den unten stehenden Link herunterladen. Old Faithful Der Inhalt des Datensatzes ist wie folgt.

1. Reihe: Letzte Eruptionsdauer (Minuten)
1. Reihe: Wartezeit bis zum nächsten Ausbruch (Minuten)

Dieses Mal verwenden wir nur die letzte Eruptionsdauer (erste Reihe) als Probe. Aus dieser Probe wurde die folgende bimodale Verteilung erhalten.

Das Gefühl, die Verteilung zu sehen Es scheint, dass eine einfache Gaußsche Verteilung nicht modelliert werden kann. Lassen Sie uns nun eine gemischte Gauß-Verteilung modellieren, die zwei Gauß-Verteilungen kombiniert. Bei der Modellierung müssen die Parameter Mittelwert $ \ mu_1, \ mu_2 $, Varianz $ \ sigma_1 ^ 2, \ sigma_2 ^ 2 $ und Mischwahrscheinlichkeit $ \ pi $ bestimmt werden. Unter der Annahme, dass die Gaußsche Verteilung $ \ phi (x | \ mu, \ sigma ^ 2) $ ist, ist die gemischte Gaußsche Verteilung durch die folgende Gleichung gegeben.

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

EM-Algorithmus

Der Mittelwert und die Varianz der Gaußschen Verteilung können durch analytische Lösung der Maximierung der Wahrscheinlichkeitsfunktion (PRML) ermittelt werden. Siehe PRML /) §2.3.4). Es ist jedoch schwierig, die Maximierung der Wahrscheinlichkeitsfunktion der gemischten Gaußschen Verteilung analytisch zu lösen, so dass die Maximierung unter Verwendung des EM-Algorithmus durchgeführt wird, der eine der Optimierungsmethoden ist. Der EM-Algorithmus ist ein Algorithmus, der zwei Schritte wiederholt, den E-Schritt und den M-Schritt.

E (Erwartung) Schritt: Berechnen Sie die Belastungsrate.
M-Schritt (Maximierung): Berechnen Sie den gewichteten Mittelwert und die Varianz.

Ausführliche Informationen zum Algorithmus finden Sie in §8.5 von Die Elemente des statistischen Lernens. Die englische Version des PDF kann kostenlos heruntergeladen werden.

Quellcode

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

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

#Durchschnitt m,Gaußsche Varianzverteilung v
def gaussian(x, m, v):
    p = math.exp(- pow(x - m, 2) / (2 * v)) / math.sqrt(2 * math.pi * v)
    return p

#E Schritt
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


#M Schritt
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


#Log Likelihood-Funktion
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

#Funktion zum Zeichnen des Histogramms
def draw_hist(xs, bins):
    plt.hist(xs, bins=bins, normed=True, alpha=0.5)

def main():
    #Lesen Sie die erste Spalte des Datensatzes
    fp = open("faithful.txt")
    data = []
    for row in fp:
        data.append(float((row.split()[0])))
    fp.close()
    #mu, vs,Stellen Sie den Anfangswert von p ein
    p = 0.5
    ms = [random.choice(data), random.choice(data)]
    vs = [np.var(data), np.var(data)]
    T = 50  #Anzahl der Iterationen
    ls = []  #Speichern Sie das Berechnungsergebnis der logarithmischen Wahrscheinlichkeitsfunktion
    #EM-Algorithmus
    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))
    #Ergebnisplot
    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()

Ausführungsergebnis

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

Die obige Abbildung zeigt die Ergebnisse der Modellierung einer Stichprobe aus einem Datensatz mit einer gemischten Gaußschen Verteilung. Die folgende Abbildung zeigt, wie sich die logarithmische Wahrscheinlichkeit erhöht, wenn der EM-Algorithmus wiederholt wird.

Andere

Anfangswert im EM-Algorithmus

Der Mittelwert war eine zufällig ausgewählte Stichprobe, die Varianz war die Varianz der Stichprobe und die Mischwahrscheinlichkeit betrug 0,5 als Anfangswert. Es gibt verschiedene Möglichkeiten, den Anfangswert auszuwählen, aber es scheint, dass dies allein zu einem Papier (?) Führt.

Über die Multi-Peak-Verteilung

Da diese Stichprobe eine bimodale Verteilung aufweist, haben wir zwei Gaußsche Verteilungen kombiniert. Natürlich können Sie drei oder mehr Gaußsche Verteilungen kombinieren, aber die Anzahl der Parameter, die Sie bestimmen müssen, nimmt zu.

Informationen zu den Abmessungen des Datensatzes

Diesmal haben wir eine eindimensionale Stichprobe verwendet, aber Sie können den EM-Algorithmus auch für multivariate Gaußsche Verteilungen verwenden. Die multivariate Gaußsche Verteilung muss signifikant mehr Parameter bestimmen als die eindimensionale Gaußsche Verteilung (Mittelwert, Kovarianzmatrix).

[PYTHON] Versuchen Sie, eine multimodale Verteilung mithilfe des EM-Algorithmus zu modellieren

Überblick

Zunächst aus der Einzelpeak-Typverteilung

Höchstwahrscheinlich Schätzung

Verfahren

Quellcode

Ausführungsergebnis

Bimodale Verteilung

Datensatz

EM-Algorithmus

Quellcode

Ausführungsergebnis

Andere

Anfangswert im EM-Algorithmus

Über die Multi-Peak-Verteilung

Informationen zu den Abmessungen des Datensatzes