👩🏽‍✈️ 🎫 😡 Wir zerlegen den EM-Algorithmus in kleine Steine 🛍️ 👩‍👧 👩🏼‍🤝‍👩🏻

In diesem Artikel werden wir, wie Sie wahrscheinlich bereits vermutet haben, über den Geräte-EM-Algorithmus sprechen. Erstens kann der Artikel für diejenigen von Interesse sein, die sich langsam der Community der Datasainisten anschließen. Das in diesem Artikel vorgestellte Material ist für diejenigen nützlicher, die kürzlich den dritten Kurs "Suche nach Datenstrukturen" im Rahmen der Spezialisierung "Maschinelles Lernen und Datenanalyse" von MIPT und Yandex belegt haben.

Das in diesem Artikel vorgestellte Material ist gewissermaßen eine Ergänzung zur ersten Trainingswoche im oben genannten Kurs, mit der Sie einige wichtige Fragen zum Funktionsprinzip des EM-Algorithmus beantworten können. Für ein besseres Verständnis des Materials ist es für unseren geschätzten Leser ratsam, Operationen mit Matrizen (Matrixmultiplikation, Finden der Determinante einer Matrix und einer inversen Matrix) ausführen zu können, die Grundlagen der Wahrscheinlichkeitstheorie und der Matstat zu verstehen und natürlich zumindest ein grundlegendes Verständnis der grundlegenden Clustering-Algorithmen zu haben und zu verstehen, was Clustering findet beim maschinellen Lernen statt. Obwohl Sie natürlich ohne dieses Wissen den Artikel lesen können, etwas und ja, es wird sicherlich klar sein :)

Nach alter Tradition wird der Artikel auch keine tiefgreifenden theoretischen Untersuchungen enthalten, sondern mit einfachen und verständlichen Beispielen gefüllt sein. Jedes nachfolgende Beispiel erklärt die Funktionsweise des EM-Algorithmus etwas tiefer als das vorherige, was uns letztendlich direkt zur Analyse des Algorithmus selbst führt. Für jedes Beispiel wird Code geschrieben. Der gesamte Code ist in Python 2.7 geschrieben, und dafür entschuldige ich mich im Voraus. Es stellte sich heraus, dass ich jetzt diese bestimmte Version verwende, aber nach dem Wechsel zu Python 3 werde ich versuchen, den Code im Artikel zu ändern.

Machen wir uns mit dem Umriss des Artikels vertraut: 1) Betrachten wir allgemein, wie der EM-Algorithmus funktioniert.

2) Wir wiederholen die Hauptpunkte aus dem Bayes-Theorem (Formel):

P_{x_{i}} (B_{j}) = \frac{P (B_{j}) \cdot P_{B_{j}} (x_{i})}{P (x_{i})}

$P_{x_i}(B_j) = \frac{P(B_j) \centerdot P_{B_j}(x_i)}{P(x_i)}$

3) Betrachten Sie das erste Beispiel des Bayes'schen Theorems, in dem alle Elemente der Formel (oder vielmehr die Wahrscheinlichkeit) enthalten sind

P (B_{j}), P_{B_{j}} (x_{i})

$P(B_j), P_{B_j}(x_i)$ ) für das Gleichheitszeichen rechts sind bekannt. In diesem Fall müssen wir nur richtig verstehen, welche Zahlen ersetzt werden sollen, und dann grundlegende Operationen ausführen.

4) Betrachten Sie das zweite Beispiel zum Bayes-Theorem. Um das Problem zu lösen, müssen wir zusätzlich die Wahrscheinlichkeitsdichte eines bestimmten Ereignisses berechnen (

x_{i}

$x_i$ ) vorbehaltlich der Hypothese (

B_{j}

$B_j$ ) -

ρ_{B_{j}} (x_{i})

$\rho_{B_j}(x_i)$ . Für Berechnungen erhalten wir Parameter einer Zufallsvariablen, nämlich die mathematische Erwartung -

μ

$\mu$ und Standardabweichung -

σ

$\sigma$ . Wahrscheinlichkeitsdichteberechnungen werden nach folgender Formel durchgeführt:

ρ_{B_{j}} (x_{i}) = \frac{1}{σ_{j} \cdot \sqrt{2 π}} \cdot \exp (- \frac{(x_{i} - μ_{j})^{2}}{2 σ_{j}^{2}})

$\rho_{B_j}(x_i) = \frac{1}{\sigma_j \centerdot \sqrt{2 \pi}} \centerdot \exp(- \frac{(x_i-\mu_j)^2}{2 \sigma_j^2})$

Bei Verwendung der obigen Formel wird davon ausgegangen, dass die Zufallsvariable eine eindimensionale Dimension hat und normal verteilt ist (mit anderen Worten, die Gaußsche oder Gaußsche Laplace-Verteilung ist eine Wahrscheinlichkeitsverteilung).

5) Betrachten Sie das dritte Beispiel, das sich vom vorherigen nur dadurch unterscheidet, dass es den Fall einer mehrdimensionalen (in unserer zweidimensionalen Version) Normalverteilung einer Zufallsvariablen berücksichtigt. In diesem Fall wird die Dichteberechnung nach folgender Formel durchgeführt:

ρ_{B_{j}} (x_{i}) = \frac{1}{(2 π)^{k / 2} \cdot ∣ Σ_{j} ∣^{1 / 2}} \cdot \exp (- \frac{1}{2} \cdot (x_{i} - μ_{j})^{T} Σ_{j}^{- 1} (x_{i} - μ_{j}))

$\rho_{B_j}(x_i) = \frac{1}{(2\pi)^{k/2} \centerdot \mid \Sigma_j\mid ^{1/2}} \centerdot \exp(- \frac{1}{2} \centerdot (x_i-\mu_j)^T \Sigma_j^{-1} (x_i-\mu_j))$

6) Schließlich modifizieren wir das dritte Beispiel so, dass die Funktionsweise des EM-Algorithmus klar demonstriert wird.

7) Abschließend vergleichen wir die Qualität des von uns abgeleiteten Algorithmus mit der Qualität des EM-Algorithmus, der in die sklearn-Bibliothek (Klasse sklearn.mixture.GaussianMixture) integriert ist.

Wenn Sie anstelle der in der Gliederung des Artikels angegebenen Formeln japanische Zeichen gesehen haben, haben Sie bitte keine Angst - Nach dem Satz von Bayes (dem zweiten Absatz der Gliederung des Papiers) werden wir solche Beispiele analysieren, die Sie wahrscheinlich zumindest intuitiv verstehen werden. Also, lasst uns gehen!

EM-Algorithmus im Allgemeinen

In dem Kurs, auf dessen Grundlage der Artikel geschrieben wurde, wird der EM-Algorithmus als eine der Clustering-Methoden vorgestellt. Mit anderen Worten, dies ist eine Methode des maschinellen Lernens ohne Lehrer, wenn wir die wahren Antworten nicht im Voraus kennen. In unserem Artikel werden wir den Algorithmus auch als Teil einer der Clustering-Methoden einer Zufallsvariablen mit einer Normalverteilung betrachten. Es sollte jedoch beachtet werden, dass der Algorithmus eine breitere Anwendung hat. Der EM-Algorithmus (englische Erwartungsmaximierung) ist eine allgemeine Methode zum Ermitteln von Schätzungen der Wahrscheinlichkeitsfunktion in Modellen mit versteckten Variablen, mit der Sie aus einer Mischung von Verteilungen (ungefähre) komplexe Wahrscheinlichkeitsverteilungen erstellen können.

Der EM-Algorithmus bei Clusterproblemen wird als iterativer Algorithmus verwendet, der bei jeder Iteration zwei Schritte ausführt:

E-Schritt. Im ersten E-Schritt wählen wir zum Beispiel zufällig versteckte Variablen aus, in unserem Fall handelt es sich um eine mathematische Erwartung -

μ

$\mu$ und Standardabweichung -

σ

$\sigma$ . Anhand der ausgewählten Variablen berechnen wir die Wahrscheinlichkeit, jedes Objekt einem bestimmten Cluster zuzuweisen. Nachfolgende E-Schritte verwenden die in den M-Schritten definierten versteckten Variablen.

M-Schritt . Im M-Schritt berechnen wir gemäß den Wahrscheinlichkeiten der Zuordnung jedes Objekts zu einem bestimmten Cluster, der im E-Schritt erhalten wurde, versteckte Variablen neu

μ

$\mu$ und

σ

$\sigma$

Iterationen werden wiederholt, bis Konvergenz auftritt. Wir können davon ausgehen, dass Konvergenz aufgetreten ist, wenn sich die Werte der verborgenen Variablen beispielsweise innerhalb einer bestimmten Konstante nicht wesentlich ändern. Im letzten Beispiel, das im Artikel behandelt wird, werden keine Konstanten festgelegt und dementsprechend werden bei jeder Iteration keine Änderungen der Werte versteckter Variablen ermittelt. Machen wir es einfacher - wir beschränken unseren Algorithmus auf eine feste Anzahl von Schritten.

Der M-Schritt ist einfach zu verstehen, und wir werden dies im letzten Beispiel sehen. Wir werden uns im Artikel auf den E-Schritt konzentrieren.

Der E-Schritt basiert auf dem Bayes'schen Theorem, auf dessen Grundlage wir die Wahrscheinlichkeit bestimmen, mit der jedes Objekt dem einen oder anderen Cluster zugeordnet wird. Erinnern wir uns, worum es in diesem Satz geht.

Erinnern Sie sich an die Bayes-Formel

P_{x_{i}} (B_{j}) = \frac{P (B_{j}) \cdot P_{B_{j}} (x_{i})}{P (x_{i})}

$P_{x_i}(B_j) = \frac{P(B_j) \centerdot P_{B_j}(x_i)}{P(x_i)}$

P_{x_{i}} (B_{j})

$P_{x_i}(B_j)$ - die Wahrscheinlichkeit, dass mit dem Ereignis

x_{i}

$x_i$ Es gab eine Hypothese

B_{j}

$B_j$

P (B_{j})

$P(B_j)$ - Wahrscheinlichkeit einer Hypothese

B_{j}

$B_j$

P_{B_{j}} (x_{i})

$P_{B_j}(x_i)$ - Eintrittswahrscheinlichkeit eines Ereignisses

x_{i}

$x_i$ vorbehaltlich der Hypothese

B_{j}

$B_j$

P (x_{i})

$P(x_i)$ - Eintrittswahrscheinlichkeit eines Ereignisses

x_{i}

$x_i$ vorbehaltlich der Erfüllung jeder Hypothese (berechnet nach der Formel der Gesamtwahrscheinlichkeit des Ereignisses)

Vorausgesetzt, dass das Ereignis

x_{i}

$x_i$ bereits passiert, die Wahrscheinlichkeit, die Hypothese zu erfüllen

B_{j}

$B_j$ nach der obigen Formel neu bewertet. Wir können das sagen

P (B_{j})

$P(B_j)$ Dies ist die a priori Wahrscheinlichkeit der Hypothese (vor dem Test) und

P_{x_{i}} (B_{j})

$P_{x_i}(B_j)$ Wird die hintere Wahrscheinlichkeit derselben Hypothese nach dem Ereignis bewertet?

x_{i}

$x_i$ das heißt, angesichts der Tatsache, dass das Ereignis

x_{i}

$x_i$ zuverlässig passiert.

Beispiel Eins nach dem Bayes-Theorem

Stellen Sie sich vor, wir hätten Teile erhalten, die auf zwei verschiedenen Maschinen hergestellt wurden. Die folgenden Merkmale der im Lager erhaltenen Produkte sind bekannt:

Insgesamt Teile - 10.000 Stück (

N = 10000

$N = 10000$ ), von denen Teile auf der ersten Maschine hergestellt wurden - 6000 Stk. (

N 1 = 6000

$N1 = 6000$ ), auf der zweiten - 4000 Stk. (

N 2 = 4000

$N2 = 4000$ )

Der Anteil der auf der ersten Maschine hergestellten Standardprodukte (d. H. Nicht fehlerhaften Produkte) beträgt 0,9, der Anteil der auf der zweiten Maschine hergestellten Standardprodukte beträgt 0,8.

Wir haben zufällig ein Teil aus den erhaltenen Teilen extrahiert, und es stellte sich als Standard heraus.

Wir müssen die Wahrscheinlichkeit bestimmen, dass:

1) das Teil auf der 1. Maschine hergestellt wird;

2) Der Artikel wurde auf der 2. Maschine hergestellt.

Entscheidung

1) Ereignis

x

$x$ In unserem Beispiel wird ein Standardprodukt extrahiert.

2) Lassen Sie uns nun die Hypothesen bestimmen

B_{j}

$B_j$ . Wir haben zwei Teile, was zwei Hypothesen bedeutet:

Hypothese

B_{1}

$B_1$ : Ein Produkt, das versehentlich entfernt wurde, wurde auf Maschine Nr. 1 hergestellt.

Hypothese

B_{2}

$B_2$ : Ein Produkt, das versehentlich entfernt wurde, wurde auf Maschine Nr. 2 hergestellt.

Dementsprechend ist die Wahrscheinlichkeit, ein auf der ersten Maschine hergestelltes Produkt zu extrahieren, gleich

P (B_{1})

$P(B_1)$ macht wieder gut

N 1 / N = 6000 / 10000 = 0.6

$N1/N = 6000/10000 = 0.6$ .

Die Wahrscheinlichkeit, ein auf einer zweiten Maschine hergestelltes Produkt zu extrahieren, beträgt

P (B_{2})

$P(B_2)$ macht wieder gut

N 2 / N = 0.4

$N2/N = 0.4$ . Und wir kümmern uns nicht um das Standardprodukt oder das defekte, es ist wichtig, aus welcher Charge es stammt.

3) Die Wahrscheinlichkeit, ein Standardprodukt aus auf Maschine Nr. 1 hergestellten Produkten zu extrahieren (dh nach der ersten Hypothese), entspricht dem Anteil der auf Maschine Nr. 1 hergestellten Standardprodukte und beträgt 0,9,

P_{B_{1}} (x_{i}) = 0.9

$P_{B_1}(x_i) = 0.9$ .

Die Wahrscheinlichkeit der Extraktion des Produkts ist ein Standardprodukt, das der Hypothese Nr. 2 unterliegt -

P_{B_{2}} (x_{i}) = 0.8

$P_{B_2}(x_i) = 0.8$

4) Bestimmen Sie die Wahrscheinlichkeit, ein Standardprodukt aus dem gesamten Produktsatz zu extrahieren -

P (x_{i})

$P(x_i)$ . In Übereinstimmung mit der Formel für die Gesamtwahrscheinlichkeit eines Ereignisses

P (x_{i}) = P (B_{1}) \cdot P_{B_{1}} (x_{i}) + P (B_{2}) \cdot P_{B_{2}} (x_{i}) = 0.6 \cdot 0.9 + 0.4 \cdot 0.8 = 0.54 + 0.32 = 0.86

$P(x_i) = P(B_1) \centerdot P_{B_1}(x_i) + P(B_2) \centerdot P_{B_2}(x_i) = 0.6 \centerdot 0.9 + 0.4 \centerdot 0.8 = 0.54 + 0.32 = 0.86$ . Hier verstehen wir das dann

0.14

$0.14$ - Dies ist die Wahrscheinlichkeit, das defekte Teil aus allen im Lager eingetroffenen Teilen zu extrahieren, und insgesamt ergibt sich 1

(0.86 + 0.14 = 1)

$(0.86 + 0.14 = 1)$ .

5) Jetzt haben wir alle Daten, um das Problem zu lösen.

5.1) Bestimmen Sie die Wahrscheinlichkeit, dass ein zufällig extrahiertes Standardteil auf einer Maschine hergestellt wird

№ 1

$№1$ ::

P_{x_{i}} (B_{1}) = \frac{P (B_{1}) \cdot P_{B_{1}} (x_{i})}{P (x_{i})} = \frac{0.6 \cdot 0.9}{0.86} \approx 0.63

$P_{x_i}(B_1) = \frac{P(B_1) \centerdot P_{B_1}(x_i)}{P(x_i)} = \frac{0.6 \centerdot 0.9}{0.86} \approx 0.63$

5.2) Wir bestimmen die Wahrscheinlichkeit, dass ein zufällig extrahiertes Standardteil auf einer Maschine hergestellt wird

№ 2

$№2$ ::

P_{x_{i}} (B_{2}) = \frac{P (B_{2}) \cdot P_{B_{2}} (x_{i})}{P (x_{i})} = \frac{0.4 \cdot 0.8}{0.86} \approx 0.37

$P_{x_i}(B_2) = \frac{P(B_2) \centerdot P_{B_2}(x_i)}{P(x_i)} = \frac{0.4 \centerdot 0.8}{0.86} \approx 0.37$

Daher haben wir die Hypothesen neu bewertet

B_{1}

$B_1$ und

B_{2}

$B_2$ . Nach einer Neubewertung bilden Hypothesen auch eine vollständige Gruppe von Ereignissen:

0.63 + 0.37 = 1

$0.63 + 0.37 = 1$ .

Nun ist es an der Zeit, mit dem zweiten Beispiel fortzufahren.

Zweites Beispiel für den Bayes-Satz unter Verwendung der Parameter der Normalverteilung einer Zufallsvariablen: mathematische Erwartung $\mu$ und Standardabweichung $\sigma$

Lassen Sie uns die Bedingungen des vorherigen Beispiels leicht ändern. Wir gehen davon aus, dass wir den Anteil der auf den Maschinen Nr. 1 und Nr. 2 hergestellten Standardprodukte nicht kennen, sondern dass wir stattdessen die durchschnittlichen Abmessungen der Durchmesser der Produkte und die Standardabweichung des Durchmessers der Produkte auf jeder Maschine kennen.

Maschine Nr. 1 produziert Teile mit einem Durchmesser von 64 mm und einer Standardabweichung von 4 mm.

Maschine Nr. 2 produziert Teile mit einem Durchmesser von 52 mm und einer Standardabweichung von 2 mm.

Eine wichtige Bedingung ist, dass der gesamte Satz von Produkten durch die Normalverteilung oder die Gaußsche Verteilung beschrieben wird.

Ansonsten sind die Bedingungen gleich, wir schreiben sie:

N 1 = 6000, μ_{1} = 64, σ_{1} = 4

$N1 = 6000, \mu_1 = 64, \sigma_1 = 4$

N 2 = 4000, μ_{2} = 52, σ_{2} = 2

$N2 = 4000, \mu_2 = 52, \sigma_2 = 2$

Wir gehen davon aus, dass es bei der Produktakzeptanz zu einem kleinen Zwischenfall gekommen ist, bei dem alle Produkte gemischt wurden.

Unsere Aufgabe ist es, die Details zu sortieren und für jedes die Wahrscheinlichkeit zu bestimmen, dass es auf Maschine Nr. 1 oder Maschine Nr. 2 hergestellt wurde. Wir gehen auch davon aus, dass das Teil auf der Maschine hergestellt wurde, deren Wahrscheinlichkeit höher ist.

Lösung

Zunächst analysieren wir den Lösungsalgorithmus. Wir können die Wahrscheinlichkeit jeder Hypothese leicht berechnen, kennen ihre Werte jedoch bereits aus einem früheren Beispiel:

P (B_{1}) = 0.6

$P(B_1)=0.6$ ,

P (B_{2}) = 0.4

$P(B_2)=0.4$ . Aber wie berechnen wir die Wahrscheinlichkeit, ein Standardprodukt für jede der verfügbaren Hypothesen einzeln zu beschlagnahmen? Alles ist einfach! Wir werden die Wahrscheinlichkeit nicht berücksichtigen, sondern die Wahrscheinlichkeitsdichte des zu entfernenden Elements mit seinem Durchmesserwert für jede Hypothese separat bestimmen .

Dazu verwenden wir die bekannte Formel:

ρ_{B_{j}} (x_{i}) = \frac{1}{σ_{j} \cdot \sqrt{2 π}} \cdot \exp (- \frac{(x_{i} - μ_{j})^{2}}{2 σ_{j}^{2}})

$\rho_{B_j}(x_i) = \frac{1}{\sigma_j \centerdot \sqrt{2 \pi}} \centerdot \exp(- \frac{(x_i-\mu_j)^2}{2 \sigma_j^2})$

x_{i}

$x_i$ - ein zufälliger Wert (in unserem Fall die tatsächliche Größe des Durchmessers des Produkts),

μ_{j}

$\mu_j$ - mathematische Erwartung von Zufallsvariablen

j

$j$ Hypothese (in unserem Fall die durchschnittliche Größe des Durchmessers des hergestellten Produkts

j

$j$ Werkzeugmaschine)

σ_{j}

$\sigma_j$ - Standardabweichung von Zufallsvariablen

j

$j$ Hypothese (in unserem Fall die Abweichung von der durchschnittlichen Größe des Durchmessers des hergestellten Produkts

j

$j$ Werkzeugmaschine).

Wenn wir also die Wahrscheinlichkeitsdichte im Zähler der Bayes-Formel geschickt durch die Wahrscheinlichkeit ersetzen, werden wir ähnliche Ersetzungen im Nenner durchführen.

Wir passen die Formeln speziell für unseren Fall an.

Die Wahrscheinlichkeit, dass das Teil mit Durchmesser Abmessungen

x_{i}

$x_i$ hergestellt auf Maschine Nr. 1 wird durch die Formel bestimmt:

P_{x_{i}} (B_{1}) = \frac{w_{1} \cdot ρ_{B_{1}} (x_{i})}{w_{1} \cdot ρ_{B_{1}} (x_{i}) + w_{2} \cdot ρ_{B_{1}} (x_{i})}

$P_{x_i}(B_1) = \frac{w_1 \centerdot \rho_{B_1}(x_i)}{w_1 \centerdot \rho_{B_1}(x_i) + w_2 \centerdot \rho_{B_1}(x_i)}$

Die Wahrscheinlichkeit, dass das Teil mit Durchmesser Abmessungen

x_{i}

$x_i$ hergestellt auf Maschine Nr. 2:

P_{x_{i}} (B_{2}) = \frac{w_{2} \cdot ρ_{B_{2}} (x_{i})}{w_{1} \cdot ρ_{B_{1}} (x_{i}) + w_{2} \cdot ρ_{B_{1}} (x_{i})}

$P_{x_i}(B_2) = \frac{w_2 \centerdot \rho_{B_2}(x_i)}{w_1 \centerdot \rho_{B_1}(x_i) + w_2 \centerdot \rho_{B_1}(x_i)}$

Beachten Sie, dass wir die Notation für die Wahrscheinlichkeit einer Hypothese durch ersetzt haben

P (B_{j})

$P(B_j)$ auf der

w_{j}

$w_j$ . Dies liegt an der entsprechenden Bezeichnung auf dem zuvor festgelegten Kurs. Weiterhin werden wir genau eine solche Notation verwenden.

Jetzt sind wir zu 100% fertig und bereit, das Problem zu lösen.

Wir schauen uns den Code an

Bibliotheken importieren

#  , 
import numpy as np
import random
import math
import matplotlib.pyplot as plt
%matplotlib inline
from sklearn.metrics import accuracy_score

Im zweiten Beispiel verwendete Funktionen

#      ,    
#      : .,  . 
def gaus_func_01(mu,sigma,x):
    return math.e**(-(x-mu)**2/(2*sigma**2)) / (sigma*(2*math.pi)**0.5)

#        
def proba_x(X, w1, w2, mu_1, mu_2, sigma_1, sigma_2):
    for i in X:
        P1_x = gaus_func_01(mu_1,sigma_1,i)
        P2_x = gaus_func_01(mu_2,sigma_2,i)
        P_x = w1*P1_x + w2*P2_x
        P_x_1 = (w1*P1_x)/P_x
        P_x_2 = (w2*P2_x)/P_x
        proba_temp = []
        proba_temp.append(P_x_1)
        proba_temp.append(P_x_2)
        proba_X.append(proba_temp)
    return proba_X

#        
def pred_x(proba_X, limit_proba):
    pred_X = []
    for x in proba_X:
        if x[0] >= limit_proba:
            pred_X.append(1)
        else:
            pred_X.append(2)
    return np.array(pred_X)

#    
def graph_01(X, pred_X, mu_1, sigma_1, mu_2, sigma_2):
    true_pred = []
    false_pred_1 = []
    false_pred_2 = []
    for i in range(X.shape[0]):
        if pred_X[i] == y[i]:
            true_pred.append([X[i], -0.025])
        else:
            if y[i] == 1:
                false_pred_1.append([X[i], -0.0075])
            else:
                false_pred_2.append([X[i], -0.015])

    false_pred_1 = np.array(false_pred_1)            
    false_pred_2 = np.array(false_pred_2)
    true_pred = np.array(true_pred)

    x_theory = np.linspace(42, 85, 20000)
    y_theory_1 = []
    for x in x_theory:
        y_theory_1.append(gaus_func_01(mu_1,sigma_1,x))
    y_theory_2 = []
    for x in x_theory:
        y_theory_2.append(gaus_func_01(mu_2,sigma_2,x))

    plt.figure(figsize=(18, 8))    
    plt.plot(
        x_theory, y_theory_1, color = 'green', lw = 2, label = 'Theoretical probability density for machine 1')
    plt.plot(
        x_theory, y_theory_2, color = 'firebrick', lw = 2, label = 'Theoretical probability density for machine 2')
    plt.hist(
        X[:N1], bins = 'auto', color='#539caf', normed = True, alpha = 0.35, label = 'machine tool products 1')
    plt.hist(
        X[N1:N], bins = 'auto', color='sandybrown', normed = True, alpha = 0.75, label = 'machine tool products 2')
    plt.plot(mu_1, 0, 'o', markersize = 11, color = 'blue', label = 'Mu 1')
    plt.plot(mu_2, 0, 'o', markersize = 11, color = 'red', label = 'Mu 2')

    plt.plot([mu_1 - sigma_1, mu_1 - sigma_1], [0,0.85*np.max(y_theory_1)],
             ':', lw = 3, color = 'blue', alpha = 0.55, label = 'Mu1 - sigma1')
    plt.plot([mu_1 + sigma_1, mu_1 + sigma_1], [0,0.85*np.max(y_theory_1)],
             ':', lw = 3, color = 'blue', alpha = 0.55, label = 'Mu1 + sigma1')
    plt.plot([mu_2 - sigma_2, mu_2 - sigma_2], [0,0.85*np.max(y_theory_2)],
             ':', lw = 3, color = 'red', alpha = 0.55, label = 'Mu2 - sigma2')
    plt.plot([mu_2 + sigma_2, mu_2 + sigma_2], [0,0.85*np.max(y_theory_2)],
             ':', lw = 3, color = 'red', alpha = 0.55, label = 'Mu2 + sigma2')

    plt.plot([mu_1 - 2 * sigma_1, mu_1 - 2 * sigma_1], [0, 0.9*0.5 * np.max(y_theory_1)],
             ':', lw = 2.5, color = 'blue', alpha = 0.35, label = 'Mu1 - 2*sigma1')
    plt.plot([mu_1 + 2 * sigma_1, mu_1 + 2 * sigma_1], [0, 0.9*0.5 * np.max(y_theory_1)],
             ':', lw = 2.5, color = 'blue', alpha = 0.35, label = 'Mu1 + 2*sigma1')
    plt.plot([mu_2 - 2 * sigma_2, mu_2 - 2 * sigma_2], [0, 0.9*0.5 * np.max(y_theory_2)],
             ':', lw = 2.5, color = 'red', alpha = 0.35, label = 'Mu2 - 2*sigma2')
    plt.plot([mu_2 + 2 * sigma_2, mu_2 + 2 * sigma_2], [0, 0.9*0.5 * np.max(y_theory_2)],
             ':', lw = 2.5, color = 'red', alpha = 0.35, label = 'Mu2 + 2*sigma2')

    plt.plot(false_pred_1[:,0], false_pred_1[:,1], 'o', markersize = 2.5, color = 'blue', alpha = 0.2, label = 'errors1')
    plt.plot(false_pred_2[:,0], false_pred_2[:,1], 'o', markersize = 2.5, color = 'red', alpha = 0.3, label = 'errors2')
    plt.plot(true_pred[:,0], true_pred[:,1], 'o', markersize = 3, color = 'green', alpha = 0.2, label = 'right answers')

    plt.xlabel('Caliber')
    plt.ylabel('Probability density')
    plt.legend()
    plt.show()

Geben Sie die Parameter an

#    
#      №1
N1 = 6000
#      №2
N2 = 4000
#      
N = N1+N2

#    №1
mu_1 = 64.
#        №1
sigma_1 = 3.5

#    №2
mu_2 = 52
#        №2
sigma_2 = 2.

Wir werden berechnen

X = np.zeros((N))
np.random.seed(seed=42)
#    ,   №1
X[:N1] = np.random.normal(loc=mu_1, scale=sigma_1, size=N1)
#  ,   №2
X[N1:N] = np.random.normal(loc=mu_2, scale=sigma_2, size=N2)

#   
y = np.zeros((N))
y[:N1] = np.array((1))
y[N1:N] = np.array((2))

#     ,    №1
w1 = float(N1)/N
#     ,    №2
w2 = float(N2)/N

#           
proba_X = []
proba_X = proba_x(X, w1, w2, mu_1, mu_2, sigma_1, sigma_2)

#   ,   ,        
limit_proba = 0.5

#     
pred_X = []
pred_X = pred_x(proba_X, limit_proba)

#    
print '   :', round(accuracy_score(y, pred_X),3)
print

print ' №1'
graph_01(X, pred_X, mu_1, sigma_1, mu_2, sigma_2)

Grafik 1

Was haben wir gerade getan? Wir haben pseudozufällig 10.000 durch die Normalverteilung beschriebene Werte erzeugt, davon 6.000 Werte mit einem Mittelwert von 64 mm und einer Standardabweichung von 4 mm, 4.000 Werte mit einem Mittelwert von 52 mm und einer Standardabweichung von 2 mm. Ferner wurde gemäß dem obigen Algorithmus für jedes Teil die Wahrscheinlichkeit seiner Herstellung auf Maschine Nr. 1 oder Nr. 2 bestimmt. Danach haben wir die Hypothese über die Herstellung des Produkts auf der einen oder anderen Maschine gewählt, je nachdem, welche Hypothese am wahrscheinlichsten ist. Schließlich haben wir die Ergebnisse unseres Algorithmus mit den wahren Antworten verglichen.

Am Ausgang haben wir einen Anteil an richtigen Antworten erhalten - 0,986. Im Allgemeinen ist dies sehr gut, da wir das Training ohne echte Antworten durchgeführt haben.

Schauen wir uns das Diagramm an. Was ist zu beachten? Sehen Sie, wo sich die Produkte befinden, die vom Algorithmus nicht korrekt identifiziert wurden.

Erstens sehen wir Algorithmusfehler nur in dem Bereich, in dem Produkte, die auf verschiedenen Maschinen hergestellt werden, den gleichen Durchmesser haben. Es sieht ziemlich logisch aus.

Zweitens sehen wir, dass der Algorithmus nur in den Objekten falsch ist, die am weitesten vom wahren Durchschnittswert der Objekte entfernt sind und gleichzeitig ziemlich nahe am falschen Zentrum liegen.

Das Interessanteste ist jedoch, dass die Probleme hauptsächlich nach dem Überschreiten ungefähr gleicher Werte beginnen

μ + - 2 σ

$\mu +- 2\sigma$ , speziell in unserem Fall an der Kreuzung

μ_{1} - 2 σ_{1}

$\mu_1 - 2\sigma_1$ und

μ_{2} + 2 σ_{2}

$\mu_2 + 2\sigma_2$

Wer möchte, kann die Parameter "verdrehen" und sehen, wie sich die Qualität des Algorithmus ändert. Dies ist nützlich für eine bessere Assimilation des Materials.

Nun, wir fahren mit dem nächsten Beispiel fort

Das dritte Beispiel. Zweidimensionaler Fall

Dieses Mal werden wir den Fall betrachten, in dem wir nicht nur Daten über den Durchmesser der Produkte haben, sondern beispielsweise auch über das Gewicht jedes Produkts. Lassen Sie Maschine Nr. 1 Teile mit einem Gewicht von 12 g und einer Standardabweichung von 1 g produzieren, Maschine Nr. 2 produziert Produkte mit einem Gewicht von 10 g und einer Standardabweichung von 0,8 mm.

Im Übrigen verwenden wir die Daten aus dem vorherigen Beispiel und schreiben sie auf.

N 1 = 6000, μ_{11} = 64, σ_{11} = 4, μ_{12} = 14, σ_{12} = 1

$N1 = 6000, \mu_{1 1} = 64, \sigma_{1 1} = 4, \mu_{1 2} = 14, \sigma_{1 2} = 1$

N 2 = 4000, μ_{21} = 52, σ_{21} = 2, μ_{22} = 10, σ_{22} = 0.9

$N2 = 4000, \mu_{2 1} = 52, \sigma_{2 1} = 2, \mu_{2 2} = 10, \sigma_{2 2} = 0.9$

Die Situation ist dieselbe - alle Details werden in einem großen Stapel gemischt, und wir müssen sie aussortieren und für jedes Teil die Wahrscheinlichkeit seiner Produktion auf Maschine Nr. 1 und Maschine Nr. 2 bestimmen.

Lösung

Im Allgemeinen ändert sich der Lösungsansatz für das vorherige Beispiel nicht, mit Ausnahme der Formel zur Bestimmung der Wahrscheinlichkeitsdichte einer Zufallsvariablen. Für den mehrdimensionalen Fall wird die folgende Formel verwendet:

p_{B_{j}} (x_{i}) = \frac{1}{(2 π)^{k / 2} \cdot ∣ Σ_{j} ∣^{1 / 2}} \cdot \exp (- \frac{1}{2} \cdot (x_{i} - μ_{j})^{T} Σ_{j}^{- 1} (x_{i} - μ_{j}))

$p_{B_j}(x_i) = \frac{1}{(2\pi)^{k/2} \centerdot \mid \Sigma_j\mid ^{1/2}} \centerdot \exp(- \frac{1}{2} \centerdot (x_i-\mu_j)^T \Sigma_j^{-1} (x_i-\mu_j))$

Σ_{j}

$\Sigma_j$ Ist eine Matrix von Kovarianzen von Zufallsvariablen

j

$j$ Hypothese (in unserem Fall die Matrix der Kovarianzen der Parameter von Produkten, die auf hergestellt wurden

j

$j$ Werkzeugmaschine)

∣ Σ_{j} ∣

$\mid \Sigma_j\mid$ Ist die Determinante der Kovarianzmatrix von Zufallsvariablen

j

$j$ Hypothese

μ_{j}

$\mu_j$ - mathematische Erwartung von Zufallsvariablen

j

$j$ Hypothese (in unserem Fall der Durchschnittswert der hergestellten Produkte

j

$j$ Werkzeugmaschine)

x_{i}

$x_i$ - Zufallsvariable (in unserem Fall die Parameter

i

$i$ -th Produkt)

Dieses Mal benötigt ein angesehener Leser Kenntnisse über Operationen mit Matrizen, um den Code zu verstehen

Im dritten Beispiel verwendete Funktionen

#      ,    
#      : .,  . 
def gaus_func_02(k, m, x, mu, sigma):
    pj_xi = []
    for j in range(k):
        det_sigma_j = np.linalg.det(sigma[j])
        factor_1 = 1 / (((2 * math.pi)**(k/2)) * ((det_sigma_j)**0.5))
        factor_2 = []
        for i in x:
#    ,            
#          
#   .    ,      
#     
            factor_2.append(math.e**float(
            -0.5 * np.matrix(i - mu[j]) * np.matrix(np.linalg.inv(sigma[j])) * np.matrix(i - mu[j]).T))
        pj_xi.append(factor_1 * np.array(factor_2))
    return np.array(pj_xi)

#     ,       №1  №2
def proba_func_02(pjxi, w, k):
    #          
    P_X = []
    for j in range(k):
        P_X.append(w[j] * pjxi[j])
    P_X = np.sum(np.array(P_X), axis = 0)
    #    ,       №1  №2
    P_J_X = []
    for j in range(k):
        P_J_X.append(w[j] * pjxi[j] / P_X)
    return np.array(P_J_X)

#        
def pred_x_02(proba_X, limit_proba):
    pred_X = []
    for x in proba_X[0]:
        if x >= limit_proba:
            pred_X.append(1)
        else:
            pred_X.append(2)
    return np.array(pred_X)

#             
def graph_02_algorithm(pred_X, mu):
    #   
    pred_X = np.array(pred_X)

    #   ,         
    answers_1 = []
    answers_2 = []

    for i in range(pred_X.shape[0]):
        if pred_X[i] == 1:
            answers_1.append(X[i])
        else:
            answers_2.append(X[i])
    
    print ' "     "'
    plt.figure(figsize=(16, 6))  
    plt.plot(
        np.array(answers_1)[:,0], np.array(answers_1)[:,1], 'o', alpha = 0.7, color='sandybrown', label = 'Produced on machine #1')
    plt.plot(
        np.array(answers_2)[:,0], np.array(answers_2)[:,1], 'o', alpha = 0.45, color = 'darkblue', label = 'Produced on machine #2')
    plt.plot(mu[0][0], mu[0][1], 'o', markersize = 16, color = 'red', label = 'Mu 1')
    plt.plot(mu[1][0], mu[1][1], 'o',  markersize = 16, color = 'slateblue', label = 'Mu 2')
    plt.xlabel('Diameter')
    plt.ylabel('Weight')
    plt.legend()
    plt.show()
    
#          
def graph_02_true(X, mu):
    print ' "  "'
    plt.figure(figsize=(16, 6))  
    plt.plot(
        X[0:N1,0], X[0:N1,1], 'o', alpha = 0.7, color='sandybrown', label = 'Produced on machine #1')
    plt.plot(
        X[N1:N,0], X[N1:N,1], 'o', alpha = 0.45, color = 'darkblue', label = 'Produced on machine #2')
    plt.plot(mu[0][0], mu[0][1], 'o', markersize = 16, color = 'red', label = 'Mu 1')
    plt.plot(mu[1][0], mu[1][1], 'o',  markersize = 16, color = 'slateblue', label = 'Mu 2')
    plt.xlabel('Diameter')
    plt.ylabel('Weight')
    plt.legend()
    plt.show()

Geben Sie die Parameter an

#  
k = 2

#      №1
N1 = 6000
#      №2
N2 = 4000
N = N1+N2

#   (  )
m = 2

#    №1
mu_1_1 = 64.
#    №1
mu_1_2 = 14.
#  
sigma_1_1 = 3.5
sigma_1_2 = 1.

#    №2
mu_2_1 = 52.
#    №2
mu_2_2 = 9.5
#  
sigma_2_1 = 2.
sigma_2_2 = 0.7

Wir werden berechnen

X = np.zeros((N, m))
np.random.seed(seed=42)
#    ,   №1
X[:N1, 0] = np.random.normal(loc=mu_1_1, scale=sigma_1_1, size=N1)
X[:N1, 1] = np.random.normal(loc=mu_1_2, scale=sigma_1_2, size=N1)
#  ,   №2
X[N1:N, 0] = np.random.normal(loc=mu_2_1, scale=sigma_2_1, size=N2)
X[N1:N, 1] = np.random.normal(loc=mu_2_2, scale=sigma_2_2, size=N2)

#    (   ,    )
y = np.zeros((N))
y[:N1] = np.array((1))
y[N1:N] = np.array((2))

#           (  )
mu  = np.array(([mu_1_1, mu_1_2], [mu_2_1, mu_2_2]))

#        (  )
sigma = np.array(([sigma_1_1, 0.],[0., sigma_1_2], [sigma_2_1, 0.],[0., sigma_2_2]))
sigma = sigma.reshape(k, m, m)

#     ,    №1  №2
w = np.array([float(1./k), float(1./k)])

#   
pj_xi = gaus_func_02(k, m, X, mu, sigma)
proba_X = proba_func_02(pj_xi, w, m)

#   ,   ,        
limit_proba = 0.5

pred_X = pred_x_02(proba_X, limit_proba)
        
#    
print '   :', round(accuracy_score(y, pred_X),3)
print

graph_02_algorithm(pred_X, mu)

graph_02_true(X, mu)

Grafik Nr. 2.1 „Verteilung der Produkte gemäß dem Algorithmus“

Grafik Nr. 2.2 „Echte Verteilung der Produkte“

In Analogie zum vorherigen Beispiel haben wir 10.000 Werte gemäß den obigen Parametern generiert

μ

$\mu$ und

σ

$\sigma$ , schrieb mehrere Funktionen für den Betrieb unseres Algorithmus und startete ihn. Der grundlegende Unterschied zwischen diesem Code und dem Code aus dem vorherigen Beispiel besteht darin, dass wir diesmal Matrixausdrücke für Berechnungen verwendet haben.

Infolgedessen zeigt unser Algorithmus den Anteil der richtigen Antworten von 0,998, was eigentlich ziemlich gut ist.

Wie wir sehen, sind die Probleme alle gleich - Fehler in den Details, die auf verschiedenen Maschinen erzeugt wurden und gleichzeitig ähnliche Abmessungen und Gewichte haben.

Der Code ist so geschrieben, dass der Leser beliebige Parameterwerte ersetzen kann

μ

$\mu$ und

σ

$\sigma$ und sehen, wie der Algorithmus funktioniert: In welchen Fällen wird sich die Qualität verschlechtern, in welchen Fällen wird sie sich verbessern. Die Hauptsache ist nicht zu vergessen, den Zeitplan zu studieren. Nun, wir kommen zu unserem letzten Punkt beim Parsen des EM-Algorithmus, tatsächlich des EM-Algorithmus selbst.

Lernen Sie den EM-Algorithmus kennen

Wir setzen unser Beispiel mit den Teilen fort, die im Lager angekommen sind. Diesmal werden wir jedoch nur wissen, dass die Produkte auf zwei verschiedenen Maschinen hergestellt wurden, es gibt 10.000 davon, jedes Teil hat einen Durchmesser und eine Größe und wir wissen nichts anderes. Aber die Aufgabe hat sich nicht geändert - wir müssen nach wie vor aus dem ganzen großen Stapel zufällig gemischter Gegenstände bestimmen, zu welcher Maschine dieser oder jener Teil gehört.

Auf den ersten Blick klingt es fast unrealistisch, aber tatsächlich haben wir ein leistungsstarkes Toolkit in der Hand: die Bayes-Formel und die Wahrscheinlichkeitsdichteformel einer Zufallsvariablen. Nehmen wir all diese Güte.

Lösung

Was werden wir tun? Wie es im EM-Algorithmus sein sollte, initialisieren wir zuerst die Parameter: Die

Wahrscheinlichkeit der Hypothese besteht darin, das auf der Maschine Nr. 1 produzierte Teil zu extrahieren -

w_{1}

$w_1$ Wir bestimmen die Wahrscheinlichkeit einer Hypothese, das auf der Maschine Nr. 2 hergestellte Teil zu extrahieren.

w_{2}

$w_2$ . Es gibt nur zwei Hypothesen, was bedeutet, dass jede von ihnen im ersten Schritt gleich 0,5 ist.

Mathematische Erwartung von Zufallsvariablen

μ

$\mu$ wie folgt definieren. Wir mischen alle Produkte mit der Zufallsfunktion, teilen die Menge gleichmäßig in zwei Teile und bestimmen für jeden Teil für jeden Parameter (Durchmesser, Gewicht) den Durchschnittswert.

Standardabweichung

σ

$\sigma$ nimm das, was von der Decke genannt wird - setze es in jeder Hinsicht gleich der Einheit. Wir schreiben im Format der Kovarianzmatrix.

Wir sind bereit, den ersten E-Schritt des Algorithmus durchzuführen. Anhand der initialisierten Parameter von Zufallsvariablen bestimmen wir die Wahrscheinlichkeit, dass jedes Teil der Maschine Nr. 1 oder der Maschine Nr. 2 zugeordnet wird.

Auf diese Weise haben wir den ersten E-Schritt gemacht.

Jetzt liegt es am M-Schritt. Hier ist alles einfach. Nachdem wir die Wahrscheinlichkeit bestimmt haben, dass jedes Teil auf einer bestimmten Maschine hergestellt wird, können wir die Wahrscheinlichkeit jeder Hypothese neu berechnen -

w_{1}

$w_1$ ,

w_{2}

$w_2$ , und auch

μ

$\mu$ und

σ

$\sigma$ .

Wir werden 15 solcher Iterationen in jeweils zwei Schritten durchführen.

Wir betrachten den Code.

Im EM-Algorithmus verwendete Funktionen

#   E-
def e_step(x, k, m, n, w, mu, sigma):
    #      i-     j-  
    pj_xi = []
    for j in range(k):
        det_sigma_j = np.linalg.det(sigma[j])
        factor_1 = 1 / (((2 * math.pi)**(k/2)) * ((det_sigma_j)**0.5))
        factor_2 = []
        for i in x:
            factor_2.append(math.e**float(
                -0.5 * np.matrix(i - mu[j]) * np.matrix(np.linalg.inv(sigma[j])) * np.matrix(i - mu[j]).T))
        pj_xi.append(factor_1 * np.array(factor_2))
    pj_xi = np.array(pj_xi)
    
    #     ,  i-    j- 
    pj_xi_w = []
    for j in range(k):
        pj_xi_w.append(pj_xi[j] * w[j])
    pj_xi_w = np.array(pj_xi_w)
    
    #     i-      
    sum_pj_xi_w = np.sum(pj_xi_w, axis = 0)
    
    #         
    proba_xi = []
    for j in range(k):
        proba_xi.append(pj_xi_w[j] / sum_pj_xi_w)
    
    return np.array(proba_xi)

#  ,    ,            ,
#       
def x_new(proba_xi):
    X1_new_ind = []
    X2_new_ind = []
    X_answers = []

    count = 0
    for x in proba_xi[0]:
        if x >= 0.5:
            X1_new_ind.append(count)
            X_answers.append(1)
        else:
            X2_new_ind.append(count)
            X_answers.append(2)
        count += 1
    #      ,    №1, №2   
    return X1_new_ind, X2_new_ind, X_answers


#   M-
def m_step(x, proba_xi,N):
    w_new = np.sum(proba_xi, axis = 1) / N
    
    #   
    mu_new = (np.array((np.matrix(proba_xi) * np.matrix(X))).T / np.sum(proba_xi, axis = 1)).T
    
    #  
    cov_new = []
    for mu in range(mu_new.shape[0]):
        X_cd = []
        X_cd_proba = []
        count = 0
        for x_i in x:
            cd = np.array(x_i - mu_new[mu])
            X_cd.append(cd)
            X_cd_proba.append(cd * proba_xi[mu][count])
            count += 1
        X_cd = np.array(X_cd)
        X_cd = X_cd.reshape(N, m)
        X_cd_proba = np.array(X_cd_proba)
        X_cd_proba = X_cd_proba.reshape(N, m)

        cov_new.append(np.matrix(X_cd.T) * np.matrix(X_cd_proba))
    cov_new = np.array((np.array(cov_new) / (np.sum(proba_xi, axis = 1)-1)))
    #             ,    
    #       ,        
    if cov_new[0][0][1] < 0:
        cov_new[0][0][1] = 0
    if cov_new[0][1][0] < 0:
        cov_new[0][1][0] = 0
    
    if cov_new[1][0][1] < 0:
        cov_new[1][0][1] = 0
    if cov_new[1][1][0] < 0:
        cov_new[1][1][0] = 0
    
    #   
    sigma_new = cov_new**0.5
    return w_new, mu_new, sigma_new

Geben Sie die Parameter an

#   
#      №1 ( 1)
N1 = 6000
#      №2 ( 2)
N2 = 4000
#       
N = N1 + N2

#  
k = 2

#    №1
mu_samples_1_1 = 64.
#    №1
mu_samples_1_2 = 14.

#    №2
mu_samples_2_1 = 52.
#    №2
mu_samples_2_2 = 9.5

#      №1
sigma_samples_1_1 = 3.5
#      №1
sigma_samples_1_2 = 1.

#      №2
sigma_samples_2_1 = 2.
#      №2
sigma_samples_2_2 = 0.7

Wir werden berechnen

X = np.zeros((N, 2))

np.random.seed(seed=42)
#    ,    №1
X[:N1, 0] = np.random.normal(loc=mu_samples_1_1, scale=sigma_samples_1_1, size=N1)
X[:N1, 1] = np.random.normal(loc=mu_samples_1_2, scale=sigma_samples_1_2, size=N1)
#    ,    №2
X[N1:N, 0] = np.random.normal(loc=mu_samples_2_1, scale=sigma_samples_2_1, size=N2)
X[N1:N, 1] = np.random.normal(loc=mu_samples_2_2, scale=sigma_samples_2_2, size=N2)

#   
m = X.shape[1]

#   
n = X.shape[0]

#        (   )
y = np.zeros((N))
y[:N1] = np.array((1))
y[N1:N] = np.array((2))

#     ,    №1  №2
w = np.array([float(1./k), float(1./k)])

np.random.seed(seed = None)
#        (   )
mu  = np.array(
    (np.mean(X[np.random.choice(n, n/k)], axis = 0), np.mean(X[np.random.choice(n, n/k)], axis = 0)))
# mu = np.array(([mu_samples_1_1, mu_samples_1_2],[mu_samples_2_1, mu_samples_2_2]))

#         (    )
sigma = np.array(([1., 0.],[0., 1.], [1., 0.],[0., 1.]))
# sigma = np.array(([sigma_samples_1_1, 0.],[0., sigma_samples_1_2], [sigma_samples_2_1, 0.],[0., sigma_samples_2_2]))
sigma = sigma.reshape(k, m, m)

#    EM-
steps = 15
#   EM-
for i in range(steps):
    proba_xi = e_step(X, k, m, n, w, mu, sigma)
    w, mu, sigma = m_step(X, proba_xi,N)
    X1_new_ind, X2_new_ind, X_answers = x_new(proba_xi)
    print ' №', i+1
    print
    print '   '
    print mu
    print
    print '   '
    print sigma
    print
    print '   '
    print round(accuracy_score(y, X_answers),3)
    
    plt.figure(figsize=(16, 6))  
    plt.plot(
        X[X1_new_ind,0], X[X1_new_ind,1], 'o', alpha = 0.7, color='sandybrown', label = 'Produced on machine #1')
    plt.plot(
        X[X2_new_ind,0], X[X2_new_ind,1], 'o', alpha = 0.45, color = 'darkblue', label = 'Produced on machine #2')
    plt.plot(mu[0][0], mu[0][1], 'o', markersize = 16, color = 'red', label = 'Mu 1')
    plt.plot(mu[1][0], mu[1][1], 'o',  markersize = 16, color = 'purple', label = 'Mu 2')
    plt.xlabel('Diameter')
    plt.ylabel('Weight')
    plt.legend()
    plt.show()

Grafik 3.1 „Verteilung der Teile bei jeder Iteration des EM-Algorithmus“

Alle Iterationen

Was ist am Ende passiert? Und es stellte sich Folgendes heraus, ohne vorläufige Informationen über die Wahrscheinlichkeit von Hypothesen für die Extraktion von Teilen, die auf Maschine Nr. 1 hergestellt wurden -

w_{1}

$w_1$ oder Maschinennummer 2

w_{2}

$w_2$ , mathematische Erwartung und Standardabweichung von Zufallsvariablen -

μ

$\mu$ und

σ

$\sigma$ haben wir eine Qualität erreicht, die mit dem vorherigen Beispiel vergleichbar ist, in dem alle oben genannten Parameter bekannt waren.

Vergleichen Sie die Parameter, die wir beim Generieren des Produktsatzes festgelegt haben, die sogenannten True-Parameter und die Parameter, die der EM-Algorithmus für uns ermittelt hat:

True-Parameter

μ_{11} = 64, μ_{12} = 14, σ_{11} = 3.5, σ_{12} = 1.

$\mu_{11} = 64, \mu_{12} = 14, \sigma_{11} = 3.5, \sigma_{12} = 1.$

μ_{21} = 52, μ_{12} = 9.5, σ_{21} = 2.0, σ_{22} = 0.7

$\mu_{21} = 52, \mu_{12} = 9.5, \sigma_{21} = 2.0, \sigma_{22} = 0.7$

Spezifische Parameter

μ_{11} \approx 64.03, μ_{12} \approx 13.99, σ_{11} \approx 3.44, σ_{12} \approx 1.23

$\mu_{11} \approx 64.03, \mu_{12} \approx 13.99, \sigma_{11} \approx 3.44, \sigma_{12} \approx 1.23$

μ_{21} \approx 52.06, μ_{12} \approx 9.54, σ_{21} \approx 1.65, σ_{22} \approx 0.77

$\mu_{21} \approx 52.06, \mu_{12} \approx 9.54, \sigma_{21} \approx 1.65, \sigma_{22} \approx 0.77$

Auf diesen Artikel geht ein Ende. Abschließend möchte ich hinzufügen, dass es für die Datenanalyse nicht sinnvoll ist, das Rad neu zu erfinden und den EM-Algorithmus selbst zu schreiben. Sie können einfach die bereitgestellte Funktion der sklearn-Bibliothek verwenden.

Mal sehen, wie die sklearn GaussianMixture-Bibliothek in unserem Beispiel funktioniert.

Modul importieren

from sklearn.mixture import GaussianMixture as GMM

Geben Sie die Parameter an

#   
#      №1 ( 1)
N1 = 6000
#      №2 ( 2)
N2 = 4000
#       
N = N1 + N2

#    №1
mu_samples_1_1 = 64.
#    №1
mu_samples_1_2 = 14.

#    №2
mu_samples_2_1 = 52.
#    №2
mu_samples_2_2 = 9.5

#      №1
sigma_samples_1_1 = 3.5
#      №1
sigma_samples_1_2 = 1.

#      №2
sigma_samples_2_1 = 2.
#      №2
sigma_samples_2_2 = 0.7


#  
k = 2

X = np.zeros((N, 2))

np.random.seed(seed=42)
#    ,    №1
X[:N1, 0] = np.random.normal(loc=mu_samples_1_1, scale=sigma_samples_1_1, size=N1)
X[:N1, 1] = np.random.normal(loc=mu_samples_1_2, scale=sigma_samples_1_2, size=N1)
#    ,    №2
X[N1:N, 0] = np.random.normal(loc=mu_samples_2_1, scale=sigma_samples_2_1, size=N2)
X[N1:N, 1] = np.random.normal(loc=mu_samples_2_2, scale=sigma_samples_2_2, size=N2)

#        (   )
y = np.zeros((N))
y[:N1] = np.array((1))
y[N1:N] = np.array((2))

Führen Sie den EM-Algorithmus aus

np.random.seed(1)
model = GMM(n_components=k, covariance_type='full')
model.fit(X)


temp_predict_X = model.predict(X)
X_answers = []
for i in range(X.shape[0]):
    if temp_predict_X[i] == 0:
        X_answers.append(1)
    else:
        X_answers.append(2)
        

print '   '
print round(accuracy_score(y, X_answers),3)

Arbeitete mit einem Knall! Keine Fehler! Wie zu erwarten ist, ist die GaussianMixture-Bibliothek des Skussian schneller und besser als unser Algorithmus. Wie wir uns erinnern, bestand der Zweck des Artikels jedoch darin, die Funktionsweise von EM-Algorithmen zu verstehen. Ich hoffe, dass das Ziel erreicht wurde :)

Informationsquellen

Literatur

1) Angewandte Statistik: Klassifikationen und Dimensionsreduktion, S.A. Ayvazyan, V.M. Buchstaber, I.S. Enyukov, L.D. Meshalkin, Moskau, Finanzen und Statistik, 1989

Internetquellen

1) Artikel 'Ein sanftes Tutorial des EM-Algorithmus und seiner Anwendung auf die Parameterschätzung für Gaußsche Gemische und Hidden-Markov-Modelle', Internationales Institut für Informatik, Jeff A. Bilmes, 1998

2) Präsentation 'Mischungsmodelle & EM-Algorithmus', NY University, David Sontag

3) Kurs 'Suche nach Datenstrukturen'

4) Formel für volle Wahrscheinlichkeits- und Bayes-Formeln

5) Funktion 'sklearn.mixture.GaussianMixture'

6) Normalverteilung, Wikipedia

7) Mehrdimensional Normalverteilung, Wikipedia

8)github.com/diefimov/MTH594_MachineLearning/blob/master/ipython/Lecture10.ipynb

9) Ein Artikel über den EM-Algorithmus

10) Ein weiterer Artikel über den EM-Algorithmus

11) Github github.com/AlexanderPetrenko83/Articles

Wir zerlegen den EM-Algorithmus in kleine Steine

EM-Algorithmus im Allgemeinen

Erinnern Sie sich an die Bayes-Formel

Beispiel Eins nach dem Bayes-Theorem

Zweites Beispiel für den Bayes-Satz unter Verwendung der Parameter der Normalverteilung einer Zufallsvariablen: mathematische Erwartung μ\muund Standardabweichung σ\sigma

Das dritte Beispiel. Zweidimensionaler Fall

Lernen Sie den EM-Algorithmus kennen

Informationsquellen

More articles:

Zweites Beispiel für den Bayes-Satz unter Verwendung der Parameter der Normalverteilung einer Zufallsvariablen: mathematische Erwartung $\mu$ und Standardabweichung $\sigma$