Nous démontons l'algorithme EM en petites briques



Dans cet article, comme vous l'avez probablement déjà deviné, nous parlerons de l'algorithme EM du périphérique. Tout d'abord, l'article peut intéresser ceux qui rejoignent lentement la communauté des dataainists. Le matériel présenté dans l'article sera plus utile à ceux qui ont récemment commencé à suivre le troisième cours "Recherche de structures de données" dans la spécialisation "Apprentissage automatique et analyse des données" du MIPT et Yandex.

Le matériel présenté dans l'article, dans un sens, est un ajout à la première semaine de formation dans le cours susmentionné, à savoir, il vous permet de répondre à certaines questions importantes concernant le principe de fonctionnement de l'algorithme EM. Pour une meilleure compréhension du matériel, il est conseillé à notre estimé lecteur de pouvoir effectuer des opérations avec des matrices (multiplication matricielle, trouver le déterminant d'une matrice et d'une matrice inverse), comprendre les bases de la théorie des probabilités et matstat, et bien sûr, avoir au moins une compréhension de base des algorithmes de clustering de base et comprendre ce que le clustering a lieu dans l'apprentissage automatique. Bien que, bien sûr, sans cette connaissance, vous pouvez lire l'article, quelque chose et oui, il sera certainement clair :)

De plus, selon l'ancienne tradition, l'article ne contiendra pas de recherche théorique approfondie, mais sera rempli d'exemples simples et compréhensibles. Chaque exemple suivant expliquera le fonctionnement de l'algorithme EM un peu plus profondément que le précédent, ce qui nous conduit finalement directement à l'analyse de l'algorithme lui-même. Pour chaque exemple, du code sera écrit. Tout le code est écrit en python 2.7, et pour cela, je m'excuse à l'avance. Il s'est avéré que maintenant j'utilise cette version particulière, mais après être passé à python 3, j'essaierai de changer le code dans l'article.



Familiarisons-nous avec les grandes lignes de l' article: 1) Considérons en termes généraux le fonctionnement de l'algorithme EM.

2) Nous répétons les principaux points du théorème de Bayes (formule):

Pxi(Bj)=P(Bj)PBj(xi)P(xi)



3) Considérons le premier exemple sur le théorème bayésien, dans lequel tous les éléments de la formule (ou plutôt la probabilité P(Bj),PBj(xi)) pour le signe égal à droite sont connus. Dans ce cas, il nous suffit de comprendre correctement les numéros où substituer, puis d'effectuer les opérations de base.

4) Considérons le deuxième exemple sur le théorème de Bayes. Pour résoudre le problème, nous devrons calculer en plus la densité de probabilité d'un certain événement (xi) sous réserve de l'hypothèse (Bj) - ρBj(xi). Pour les calculs, on nous donnera les paramètres d'une variable aléatoire, à savoir l'espérance mathématique -μ et écart type - σ. Les calculs de densité de probabilité seront effectués selon la formule suivante:

ρBj(xi)=1σj2πexp((xiμj)22σj2)


L'utilisation de la formule ci-dessus suppose que la variable aléatoire a une dimension unidimensionnelle et est distribuée normalement (en d'autres termes, la distribution gaussienne ou gaussienne-Laplace est une distribution de probabilité).

5) Considérons le troisième exemple, qui ne diffère du précédent que par le fait qu'il considère le cas d'une distribution normale multidimensionnelle (dans notre version bidimensionnelle) d'une variable aléatoire. Dans ce cas, le calcul de la densité est effectué selon la formule suivante:

ρBj(xi)=1(2π)k/2Σj1/2exp(12(xiμj)TΣj1(xiμj))



6) Enfin, nous modifions le troisième exemple de manière à démontrer clairement le fonctionnement de l'algorithme EM.

7) En conclusion, nous comparons la qualité de l'algorithme que nous avons déduit avec la qualité de l'algorithme EM, qui est intégré dans la bibliothèque sklearn (classe sklearn.mixture.GaussianMixture)

Si, au lieu des formules spécifiées dans le plan de l'article, vous avez vu des caractères japonais, alors s'il vous plaît, n'ayez pas peur - après le théorème de Bayes (le deuxième paragraphe du plan de l'article), nous analyserons de tels exemples que vous comprendrez probablement au moins intuitivement. Eh bien, allons-y!

Algorithme EM en général


Sur le cours, sur la base duquel l'article a été écrit, l'algorithme EM est présenté comme l'une des méthodes de clustering. En d'autres termes, il s'agit d'une méthode d'apprentissage automatique sans professeur, lorsque nous ne connaissons pas les vraies réponses à l'avance. Dans notre article, nous considérerons également l'algorithme comme faisant partie d'une des méthodes de clustering d'une variable aléatoire avec une distribution normale. Cependant, il convient de noter que l'algorithme a une application plus large. L'algorithme EM (anglais expectation-maximization) est une méthode générale pour trouver des estimations de la fonction de vraisemblance dans des modèles avec des variables cachées, qui à partir d'un mélange de distributions vous permet de construire des distributions de probabilité complexes (approximatives).

L'algorithme EM dans les problèmes de clustering est utilisé comme un algorithme itératif, qui effectue deux étapes à chaque itération:

E-step. À la première étape E, d'une certaine manière, par exemple, au hasard, nous sélectionnons des variables cachées, dans notre cas, ce sera une attente mathématique -μ et écart type - σ. En utilisant les variables sélectionnées, nous calculons la probabilité d'affecter chaque objet à un cluster particulier. Les étapes E suivantes utilisent les variables cachées définies dans les étapes M.

M-étape . À l'étape M, nous, conformément aux probabilités d'affecter chaque objet à un cluster particulier obtenu à l'étape E, recalculons les variables cachéesμ et σ

Les itérations sont répétées jusqu'à ce que la convergence se produise. Nous pouvons supposer que la convergence s'est produite lorsque les valeurs des variables cachées ne changent pas de manière significative, par exemple, à l'intérieur d'une constante donnée. Dans le dernier exemple, qui est considéré dans l'article, nous ne définirons pas de constantes et, par conséquent, nous ne déterminerons pas les changements dans les valeurs des variables cachées à chaque itération. Faisons-le plus simplement - nous limiterons notre algorithme à un nombre fixe d'étapes.

L'étape M est assez simple à comprendre, et nous le verrons dans le dernier exemple. Nous nous concentrerons sur l'étape E dans l'article.

L'étape E est basée sur le théorème bayésien, sur la base duquel nous déterminons la probabilité que chaque objet soit assigné à l'un ou l'autre cluster. Rappelons-nous de quoi parle ce théorème.

Rappelez-vous la formule Bayes


Pxi(Bj)=P(Bj)PBj(xi)P(xi)



Pxi(Bj)- la probabilité qu'avec l'événement xiil y avait une hypothèse Bj

P(Bj)- probabilité d'hypothèse Bj

PBj(xi)- probabilité d'occurrence d'un événement xisous réserve d'hypothèse Bj

P(xi)- probabilité d'occurrence d'un événement xisous réserve de la réalisation de chaque hypothèse (calculée par la formule de la probabilité totale de l'événement)

A condition que l'événementxi déjà eu lieu, la probabilité de réaliser l'hypothèse Bjréévalué en utilisant la formule ci-dessus. On peut dire çaP(Bj) c'est la probabilité a priori de l'hypothèse (avant le test), et Pxi(Bj)La probabilité postérieure de la même hypothèse est-elle évaluée après l'événement xic'est-à-dire, étant donné que l'événement xiest arrivé de manière fiable.

Exemple un par le théorème de Bayes


Imaginez que nous recevions des pièces produites sur deux machines différentes. Les caractéristiques suivantes des produits reçus à l'entrepôt sont connues:

Au total, pièces - 10 000 pièces (N=10000), dont pièces produites sur la première machine - 6000 pièces. (N1=6000), sur le second - 4000 pcs. (N2=4000)

La proportion de produits standard (c'est-à-dire non défectueux) fabriqués sur la première machine est de 0,9, la part des produits standard fabriqués sur la deuxième machine est de 0,8.

Nous avons extrait au hasard une partie des pièces reçues, et cela s'est avéré être standard.

Nous devons déterminer la probabilité que:

1) la pièce soit produite sur la 1ère machine;

2) l'article a été produit sur la 2e machine.

Décision

1) Événementxdans notre exemple, extraire un produit standard.

2) Maintenant, déterminons les hypothèsesBj. Nous avons deux lots de parties, ce qui signifie deux hypothèses:

HypothèseB1: un produit retiré accidentellement a été fabriqué sur la machine n ° 1.

HypothèseB2: un produit retiré accidentellement a été fabriqué sur la machine n ° 2.

En conséquence, la probabilité d'extraire un produit fabriqué sur la première machine estP(B1) compose N1/N=6000/10000=0.6.

La probabilité d'extraire un produit fabriqué sur une deuxième machine estP(B2) compose N2/N=0.4. Et nous ne nous soucions pas du produit standard ou du produit défectueux, il est important de quel lot il s'agit.

3) La probabilité d'extraire un produit standard de produits fabriqués sur la machine n ° 1 (c'est-à-dire selon la première hypothèse) correspond à la proportion de produits standard produits sur la machine n ° 1 et est de 0,9,PB1(xi)=0.9.

La probabilité d'extraire le produit est un produit standard soumis à l'hypothèse n ° 2 -PB2(xi)=0.8

4) Déterminer la probabilité d'extraire un produit standard de l'ensemble des produits -P(xi). Conformément à la formule de la probabilité totale d'un événementP(xi)=P(B1)PB1(xi)+P(B2)PB2(xi)=0.60.9+0.40.8=0.54+0.32=0.86. Ici, nous comprenons que0.14 - c'est la probabilité d'extraire la pièce défectueuse de toutes les pièces qui sont arrivées à l'entrepôt et au total il s'avère 1 (0.86+0.14=1).

5) Nous avons maintenant toutes les données pour résoudre le problème.

5.1) Déterminer la probabilité qu'une pièce standard extraite au hasard soit produite sur une machine1:

Pxi(B1)=P(B1)PB1(xi)P(xi)=0.60.90.860.63



5.2) Nous déterminons la probabilité qu'une pièce standard extraite au hasard soit produite sur une machine2:

Pxi(B2)=P(B2)PB2(xi)P(xi)=0.40.80.860.37



Ainsi, nous avons réévalué les hypothèses B1et B2. Après réévaluation, les hypothèses forment également un groupe complet d'événements:0.63+0.37=1.

Eh bien, il est maintenant temps de passer au deuxième exemple.

Deuxième exemple pour le théorème de Bayes utilisant les paramètres de la distribution normale d'une variable aléatoire: l'espérance mathématique μet écart type σ



Modifions légèrement les conditions de l'exemple précédent. Nous supposons que nous ne connaissons pas la proportion de produits standard fabriqués sur les machines n ° 1 et n ° 2, mais supposons plutôt que nous connaissons les dimensions moyennes des diamètres des produits et l'écart type du diamètre des produits sur chaque machine.

La machine n ° 1 produit des pièces de 64 mm de diamètre et un écart type de 4 mm.

La machine n ° 2 produit des pièces de 52 mm de diamètre et un écart type de 2 mm.

Une condition importante est que l'ensemble des produits soit décrit par la distribution normale ou distribution gaussienne.

Sinon, les conditions sont les mêmes, nous les écrivons:

N1=6000,μ1=64,σ1=4

N2=4000,μ2=52,σ2=2

Nous supposons que lors du processus d'acceptation du produit, il y a eu un petit incident, à la suite duquel tous les produits ont été mélangés.

Notre tâche est de trier les détails et pour chacun de déterminer la probabilité qu'elle ait été produite sur la machine n ° 1 ou la machine n ° 2. Nous supposerons également que la pièce a été produite sur la machine, dont la probabilité est plus élevée.

Solution

Tout d'abord, nous analyserons l'algorithme de solution. Nous pouvons facilement calculer la probabilité de chaque hypothèse, cependant, nous connaissons déjà leurs valeurs à partir d'un exemple passé:P(B1)=0.6, P(B2)=0.4. Mais comment calculer la probabilité de saisir un produit standard individuellement pour chacune des hypothèses disponibles? Tout est simple! Nous ne considérerons pas la probabilité, au lieu de cela, nous déterminerons la densité de probabilité de l'élément à supprimer avec sa valeur de diamètre pour chaque hypothèse séparément.

Pour ce faire, nous utilisons la formule bien connue:

ρBj(xi)=1σj2πexp((xiμj)22σj2)



xi- une valeur aléatoire (dans notre cas, la taille réelle du diamètre du produit),

μj- attente mathématique des variables aléatoires jhypothèse (dans notre cas, la taille moyenne du diamètre du produit fabriqué sur jmachine-outils)

σj- écart type des variables aléatoires jhypothèse (dans notre cas, l'écart par rapport à la taille moyenne du diamètre du produit fabriqué sur jmachine-outils).

Ainsi, nous effectuons une substitution intelligente de probabilité pour la densité de probabilité dans le numérateur de la formule de Bayes, nous effectuerons des remplacements similaires dans le dénominateur.

Nous adaptons les formules spécifiquement pour notre cas.

La probabilité que la pièce avec des dimensions de diamètrexi produite sur la machine n ° 1 est déterminée par la formule:

Pxi(B1)=w1ρB1(xi)w1ρB1(xi)+w2ρB1(xi)



La probabilité que la pièce avec des dimensions de diamètre xiproduite sur la machine n ° 2:

Pxi(B2)=w2ρB2(xi)w1ρB1(xi)+w2ρB1(xi)



Notez que nous avons remplacé la notation de la probabilité d'une hypothèse par P(Bj)sur le wj. Cela est dû à la désignation correspondante sur le parcours précédemment désigné. De plus, nous utiliserons une telle notation.

Maintenant, nous sommes à 100% complets et prêts à résoudre le problème.

Nous regardons le code

Importer des bibliothèques
#  , 
import numpy as np
import random
import math
import matplotlib.pyplot as plt
%matplotlib inline
from sklearn.metrics import accuracy_score


Fonctions utilisées dans le deuxième exemple
#      ,    
#      : .,  . 
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()


Spécifiez les paramètres
#    
#      №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.


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


Graphique # 1



Que venons-nous de faire? Nous avons généré de manière pseudo-aléatoire 10000 valeurs décrites par la distribution normale, dont 6000 valeurs avec une moyenne de 64 mm et un écart type de 4 mm, 4000 valeurs avec une moyenne de 52 mm et un écart type de 2 mm. De plus, conformément à l'algorithme ci-dessus, pour chaque pièce, la probabilité de sa production sur la machine n ° 1 ou n ° 2 a été déterminée. Après cela, nous avons choisi l'hypothèse de production du produit sur l'une ou l'autre machine, en fonction de l'hypothèse la plus probable. Enfin, nous avons comparé les résultats de notre algorithme avec les vraies réponses.

À la sortie, nous avons obtenu une part de bonnes réponses - 0,986. En général, c'est très bien, étant donné que nous avons mené la formation sans utiliser de vraies réponses.

Regardons le graphique. À quoi est-il recommandé de faire attention? Voir où se trouvent les produits qui ne sont pas correctement identifiés par l'algorithme.

Premièrement, nous ne voyons des erreurs d'algorithme que dans la zone où les produits fabriqués sur différentes machines ont le même diamètre. Cela semble assez logique.

Deuxièmement, nous voyons que l'algorithme ne se trompe que dans les objets les plus éloignés de la vraie valeur moyenne des objets et en même temps assez proches du faux centre.

Mais la chose la plus intéressante est que les problèmes commencent principalement après avoir croisé des valeurs approximativement égalesμ+2σ, en particulier dans notre cas à l'intersection μ12σ1et μ2+2σ2

Ceux qui le souhaitent peuvent "tordre" les paramètres et voir comment la qualité de l'algorithme change. Cela sera utile pour une meilleure assimilation du matériau.

Eh bien, nous passons à l'exemple suivant

Le troisième exemple. Boîtier bidimensionnel


Cette fois, nous considérerons le cas lorsque nous disposerons de données non seulement sur le diamètre des produits, mais aussi, par exemple, sur le poids de chaque produit. Laissez la machine n ° 1 produire des pièces pesant 12 g et un écart-type de 1 g, la machine n ° 2 produit des produits pesant 10 g et un écart-type de 0,8 mm.

Pour le reste, nous utiliserons les données de l'exemple précédent, nous les noterons.

N1=6000,μ11=64,σ11=4,μ12=14,σ12=1

N2=4000,μ21=52,σ21=2,μ22=10,σ22=0.9

La situation est la même - tous les détails ont été mélangés en un seul gros tas et nous devons les trier et déterminer pour chaque pièce la probabilité de sa production sur la machine n ° 1 et la machine n ° 2.

Solution

En général, l'approche de la solution de l'exemple précédent ne change pas, à l'exception de la formule pour déterminer la densité de probabilité d'une variable aléatoire. Pour le cas multidimensionnel, la formule suivante est utilisée:

pBj(xi)=1(2π)k/2Σj1/2exp(12(xiμj)TΣj1(xiμj))



ΣjEst une matrice de covariances de variables aléatoires jhypothèse (dans notre cas, la matrice des covariances des paramètres des produits fabriqués sur jmachine-outils)

ΣjEst le déterminant de la matrice de covariance des variables aléatoires jhypothèse

μj- attente mathématique des variables aléatoires jhypothèse (dans notre cas, la valeur moyenne des produits fabriqués sur jmachine-outils)

xi- variable aléatoire (dans notre cas, les paramètres i-th produit)

Cette fois, pour comprendre le code, un lecteur respecté aura besoin de connaissances sur les opérations avec les matrices

Fonctions utilisées dans le troisième exemple
#      ,    
#      : .,  . 
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()


Spécifiez les paramètres
#  
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


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


Graphique n ° 2.1 «Distribution des produits selon l'algorithme»


Graphique n ° 2.2 «Distribution réelle des produits»


Par analogie avec l'exemple précédent, nous avons généré 10 000 valeurs conformément aux paramètres ci-dessusμ et σ, a écrit plusieurs fonctions pour le fonctionnement de notre algorithme et l'a lancé. La différence fondamentale entre ce code et le code de l'exemple précédent est que cette fois nous avons utilisé des expressions matricielles pour les calculs.

En conséquence, notre algorithme montre la proportion de réponses correctes égale à 0,998, ce qui est en fait assez bon.

Les problèmes, comme nous le voyons, sont tous les mêmes - des erreurs dans les détails qui ont été produites sur différentes machines et en même temps ont des dimensions et des poids similaires.

Le code est écrit pour que le lecteur puisse remplacer toutes les valeurs de paramètresμ et σet voyez comment l'algorithme fonctionnera: dans quels cas la qualité se détériorera, dans laquelle elle s'améliorera. L'essentiel est de ne pas oublier d'étudier le planning. Eh bien, nous passons à notre point final en analysant l'algorithme EM, en fait l'algorithme EM lui-même.

Rencontrez l'algorithme EM


Nous continuons notre exemple avec les pièces arrivées à l'entrepôt. Mais cette fois, nous saurons seulement que les produits ont été fabriqués sur deux machines différentes, il y en a 10 000, chaque pièce a un diamètre et une taille, et nous ne savons rien d'autre. Mais la tâche n'a pas changé - nous, comme auparavant, à partir de la grande pile entière, des éléments mélangés au hasard, nous devrons déterminer à quelle machine telle ou telle partie appartient.

À première vue, cela semble presque irréaliste, mais en fait, nous avons entre nos mains une puissante boîte à outils: la formule de Bayes et la formule de densité de probabilité d'une variable aléatoire. Prenons toute cette bonté.

Solution

Que ferons-nous? Comme il se doit dans l'algorithme EM, on initialise d'abord les paramètres: La

probabilité de l'hypothèse est d'extraire la pièce produite sur la machine n ° 1 -w1 on détermine la probabilité d'une hypothèse pour extraire la pièce produite sur la machine n ° 2 - w2. Il n'y a que deux hypothèses, ce qui signifie que chacune d'elles dans la première étape sera égale à 0,5.

Attente mathématique des variables aléatoiresμdéfinir comme suit. Nous mélangeons tous les produits en utilisant la fonction aléatoire, divisons l'ensemble de manière égale en deux parties, pour chaque partie pour chaque paramètre (diamètre, poids), nous déterminons la valeur moyenne.

Écart-typeσprenez ce qu'on appelle du plafond - mettez-le égal à tous égards à l'unité. Nous écrivons au format de la matrice de covariance.

Nous sommes prêts à effectuer la première étape E de l'algorithme. En utilisant les paramètres initialisés des variables aléatoires, nous déterminons la probabilité que chaque pièce soit affectée à la machine n ° 1 ou à la machine n ° 2.

En fait, de cette façon, nous avons fait la première étape E.

Maintenant, c'est à l'étape M. Ici, tout est simple. Après avoir déterminé la probabilité que chaque pièce soit produite sur une machine particulière, nous pouvons recalculer la probabilité de chaque hypothèse -w1, w2, et μet σ.

Nous allons effectuer 15 itérations de ce type, en deux étapes chacune.

Nous examinons le code.
Fonctions utilisées dans l'algorithme EM
#   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


Spécifiez les paramètres
#   
#      №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


Nous calculerons
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()


Graphique 3.1 «Répartition des pièces à chaque itération de l'algorithme EM»

Toutes les itérations
















Ce qui est arrivé à la fin? Et il s'est avéré ce qui suit, n'ayant aucune information préliminaire sur la probabilité d'hypothèses pour l'extraction de pièces produites sur la machine n ° 1 -w1 ou numéro de machine 2 w2, espérance mathématique et écart-type des variables aléatoires - μet σ, nous avons atteint une qualité comparable à l'exemple précédent, dans laquelle tous les paramètres ci-dessus étaient connus.

Comparez les paramètres que nous définissons lors de la génération de l'ensemble de produits, les soi-disant vrais paramètres et les paramètres que l'algorithme EM a déterminés pour nous:

vrais paramètres

μ11=64,μ12=14,σ11=3.5,σ12=1.

μ21=52,μ12=9.5,σ21=2.0,σ22=0.7

Paramètres spécifiques

μ1164.03,μ1213.99,σ113.44,σ121.23

μ2152.06,μ129.54,σ211.65,σ220.77

Sur cet article prend fin. En conclusion, j'ajouterai qu'à des fins d'analyse de données, cela n'a pas de sens de réinventer la roue et d'écrire un algorithme EM vous-même, vous pouvez simplement utiliser la fonction de bibliothèque sklearn fournie.

Voyons comment la bibliothèque sklearn GaussianMixture fonctionne sur notre exemple.

Module d'importation
from sklearn.mixture import GaussianMixture as GMM


Spécifiez les paramètres
#   
#      №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))


Exécutez l'algorithme EM
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)


A travaillé avec un bang! Pas d'erreurs! Comme vous vous en doutez, la bibliothèque GaussianMixture du skussian est plus rapide et meilleure que notre algorithme. Cependant, le but de cet article, comme nous nous en souvenons, était de comprendre le fonctionnement des algorithmes EM. J'espère que l'objectif a été atteint :)

Sources d'information


Littérature

1) Statistiques appliquées: classifications et réduction de la dimensionnalité, S.A. Ayvazyan, V.M. Buchstaber, I.S. Enyukov, L.D. Meshalkin, Moscou, Finance and Statistics, 1989

Internet resources

1) Article 'A Gentle Tutorial of the EM Algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models', International Computer Science Institute, Jeff A. Bilmes, 1998

2) Présentation 'Modèles de mélange et algorithme EM', Université de NY, David Sontag

3) Cours 'Recherche de structures de données'

4) Formule de probabilité complète et formules de Bayes

5) Fonction 'sklearn.mixture.GaussianMixture'

6) Distribution normale, Wikipedia

7) Multidimensionnel distribution normale, Wikipedia

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

9) Un article sur l'algorithme EM

10) Un autre article sur l'algorithme EM

11) Github github.com/AlexanderPetrenko83/Articles

All Articles