Desmontamos o algoritmo EM em pequenos tijolos



Neste artigo, como você provavelmente já adivinhou, falaremos sobre o algoritmo EM do dispositivo. Primeiro de tudo, o artigo pode ser de interesse para aqueles que estão se juntando lentamente à comunidade de dataainists. O material apresentado no artigo será mais útil para quem começou recentemente o terceiro curso "Pesquisando estruturas de dados" na especialização "Machine Learning e Data Analysis" do MIPT e Yandex.

O material apresentado no artigo, em certo sentido, é um complemento à primeira semana de treinamento no curso mencionado, ou seja, permite responder a algumas questões importantes sobre o princípio de operação do algoritmo EM. Para uma melhor compreensão do material, é aconselhável que o nosso estimado leitor seja capaz de executar operações com matrizes (multiplicação de matrizes, descobrindo o determinante de uma matriz e uma matriz inversa), entender os conceitos básicos da teoria das probabilidades e do matstat e, é claro, ter pelo menos um entendimento básico dos algoritmos básicos de agrupamento e entender o que o armazenamento em cluster ocorre no aprendizado de máquina. Embora, é claro, sem esse conhecimento, você possa ler o artigo, algo e sim certamente ficará claro :)

Além disso, de acordo com a antiga tradição, o artigo não conterá uma pesquisa teórica profunda, mas será preenchido com exemplos simples e compreensíveis. Cada exemplo subseqüente explicará a operação do algoritmo EM um pouco mais profundo do que o anterior, o que nos leva, em última análise, à análise do próprio algoritmo. Para cada exemplo, o código será gravado. Todo o código é escrito em python 2.7, e por isso peço desculpas antecipadamente. Acontece que agora estou usando essa versão específica, mas depois de mudar para o python 3, tentarei alterar o código no artigo.



Vamos nos familiarizar com o esboço do artigo: 1) Vamos considerar em termos gerais como o algoritmo EM funciona.

2) Repetimos os pontos principais do teorema de Bayes (fórmula):

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



3) Considere o primeiro exemplo do teorema bayesiano, no qual todos os elementos da fórmula (ou melhor, probabilidade P(Bj),PBj(xi)) para o sinal de igual à direita. Nesse caso, precisamos apenas entender corretamente quais números substituir e executar operações básicas.

4) Considere o segundo exemplo no teorema de Bayes. Para resolver o problema, precisaremos calcular adicionalmente a densidade de probabilidade de um determinado evento (xi) sujeito à hipótese (Bj) - ρBj(xi). Para os cálculos, receberemos parâmetros de uma variável aleatória, a saber, expectativa matemática -μ e desvio padrão - σ. Os cálculos de densidade de probabilidade serão realizados de acordo com a seguinte fórmula:

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


O uso da fórmula acima pressupõe que a variável aleatória tenha uma dimensão unidimensional e seja distribuída normalmente (em outras palavras, a distribuição Gaussiana ou Gaussiana-Laplace é uma distribuição de probabilidade).

5) Considere o terceiro exemplo, que difere do anterior apenas por considerar o caso de uma distribuição normal multidimensional (em nossa versão bidimensional) de uma variável aleatória. Nesse caso, o cálculo da densidade é realizado de acordo com a seguinte fórmula:

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



6) Finalmente, modificamos o terceiro exemplo de maneira a demonstrar claramente a operação do algoritmo EM.

7) Em conclusão, comparamos a qualidade do algoritmo que deduzimos com a qualidade do algoritmo EM, que está embutido na biblioteca sklearn (classe sklearn.mixture.GaussianMixture)

Se, em vez das fórmulas especificadas no esboço do artigo, você viu caracteres japoneses, não tenha medo - após o teorema de Bayes (o segundo parágrafo do esboço do artigo), analisaremos esses exemplos que você provavelmente entenderá pelo menos intuitivamente. Bem, vamos!

Algoritmo EM em geral


No curso, com base no qual o artigo foi escrito, o algoritmo EM é apresentado como um dos métodos de agrupamento. Em outras palavras, esse é um método de aprendizado de máquina sem professor, quando não sabemos as respostas verdadeiras com antecedência. Em nosso artigo, também consideraremos o algoritmo como parte de um dos métodos de agrupamento de uma variável aleatória com uma distribuição normal. No entanto, deve-se notar que o algoritmo tem uma aplicação mais ampla. O algoritmo EM (maximização de expectativa em inglês) é um método geral para encontrar estimativas da função de verossimilhança em modelos com variáveis ​​ocultas, o que permite construir (aproximadas) distribuições de probabilidade complexas a partir de uma mistura de distribuições.

O algoritmo EM nos problemas de clustering é usado como um algoritmo iterativo, que executa duas etapas em cada iteração:

E-step. Na primeira etapa do E, de alguma forma, por exemplo, aleatoriamente, selecionamos variáveis ​​ocultas; no nosso caso, será uma expectativa matemática -μ e desvio padrão - σ. Usando as variáveis ​​selecionadas, calculamos a probabilidade de atribuir cada objeto a um cluster específico. As etapas E subsequentes usam as variáveis ​​ocultas definidas nas etapas M.

H-passo . Na etapa M, nós, de acordo com as probabilidades de atribuir cada objeto a um cluster específico obtido na etapa E, recalculamos as variáveis ​​ocultasμ e σ

As iterações são repetidas até a convergência ocorrer. Podemos assumir que a convergência ocorreu quando os valores das variáveis ​​ocultas não mudam significativamente, por exemplo, dentro de uma determinada constante. No último exemplo, que é considerado no artigo, não definiremos constantes e, portanto, não determinaremos alterações nos valores de variáveis ​​ocultas a cada iteração. Vamos simplificar - limitaremos nosso algoritmo a um número fixo de etapas.

O passo M é simples o suficiente para entender, e veremos isso no último exemplo. Vamos nos concentrar na etapa E deste artigo.

A etapa E é baseada no teorema bayesiano, com base no qual determinamos a probabilidade de cada objeto ser atribuído a um ou outro aglomerado. Vamos lembrar do que trata esse teorema.

Lembre-se da fórmula de Bayes


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



Onde Pxi(Bj)- a probabilidade de que com o evento xihavia uma hipótese Bj

P(Bj)- probabilidade de hipótese Bj

PBj(xi)- probabilidade de ocorrência de um evento xisujeito a hipótese Bj

P(xi)- probabilidade de ocorrência de um evento xisujeito ao cumprimento de cada hipótese (calculado pela fórmula da probabilidade total do evento)

Desde que o eventoxi já aconteceu, a probabilidade de cumprir a hipótese Bjreavaliado usando a fórmula acima. Nós podemos dizer queP(Bj) essa é a probabilidade a priori da hipótese (antes do teste) e Pxi(Bj)A probabilidade posterior da mesma hipótese é avaliada após o evento xiisto é, dado que o evento xiaconteceu de forma confiável.

Exemplo Um Pelo Teorema de Bayes


Imagine que recebemos peças que foram produzidas em duas máquinas diferentes. As seguintes características dos produtos recebidos no armazém são conhecidas:

No total, peças - 10.000 peças (N=10000), das quais peças produzidas na primeira máquina - 6000 peças. (N1=6000), no segundo - 4000 pcs. (N2=4000)

A proporção de produtos padrão (ou seja, não defeituosos) fabricados na primeira máquina é de 0,9, a parcela de produtos padrão fabricados na segunda máquina é de 0,8.

Extraímos aleatoriamente uma parte das partes recebidas, e acabou sendo padrão.

Precisamos determinar a probabilidade de que:

1) a peça seja produzida na 1ª máquina;

2) o item foi produzido na 2ª máquina.

Decisão

1) Eventoxem nosso exemplo, extraindo um produto padrão.

2) Agora, vamos determinar as hipótesesBj. Temos duas partes, o que significa duas hipóteses:

HipóteseB1: um produto que foi removido acidentalmente foi produzido na máquina nº 1.

HipóteseB2: um produto que foi removido acidentalmente foi produzido na máquina nº 2.

Consequentemente, a probabilidade de extrair um produto fabricado na primeira máquina éP(B1) inventa N1/N=6000/10000=0.6.

A probabilidade de extrair um produto fabricado em uma segunda máquina éP(B2) inventa N2/N=0.4. E não nos importamos com o produto padrão ou com defeito, é importante de qual lote é.

3) A probabilidade de extrair um produto padrão de produtos fabricados na máquina nº 1 (ou seja, sob a primeira hipótese) corresponde à proporção de produtos padrão produzidos na máquina nº 1 e é de 0,9,PB1(xi)=0.9.

A probabilidade de extrair o produto é um produto padrão sujeito à hipótese nº 2 -PB2(xi)=0.8

4) Determine a probabilidade de extrair um produto padrão de todo o conjunto de produtos -P(xi). De acordo com a fórmula para a probabilidade total de um eventoP(xi)=P(B1)PB1(xi)+P(B2)PB2(xi)=0.60.9+0.40.8=0.54+0.32=0.86. Aqui entendemos que então0.14 - esta é a probabilidade de extrair a peça defeituosa de todas as peças que chegaram ao armazém e, no total, resulta 1 (0.86+0.14=1).

5) Agora temos todos os dados para resolver o problema.

5.1) Determinar a probabilidade de uma peça padrão extraída aleatoriamente ser produzida em uma máquina1:

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



5.2) Determinamos a probabilidade de que uma peça padrão extraída aleatoriamente seja produzida em uma máquina2:

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



Assim, reavaliamos as hipóteses B1e B2. Após a reavaliação, as hipóteses também formam um grupo completo de eventos:0.63+0.37=1.

Bem, agora é a hora de passar para o segundo exemplo.

Segundo exemplo para o teorema de Bayes usando os parâmetros da distribuição normal de uma variável aleatória: expectativa matemática μe desvio padrão σ



Vamos modificar um pouco as condições do exemplo anterior. Assumimos que não sabemos a proporção de produtos padrão fabricados nas máquinas nº 1 e nº 2, mas supomos que, em vez disso, sabemos as dimensões médias dos diâmetros dos produtos e o desvio padrão do diâmetro dos produtos em cada máquina.

A máquina nº 1 produz peças com 64 mm de diâmetro e um desvio padrão de 4 mm.

A máquina nº 2 produz peças com 52 mm de diâmetro e um desvio padrão de 2 mm.

Uma condição importante é que todo o conjunto de produtos seja descrito pela distribuição normal ou distribuição gaussiana.

Caso contrário, as condições são as mesmas, nós as escrevemos:

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

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

Assumimos que no processo de aceitação do produto houve um pequeno incidente, como resultado do qual todos os produtos foram misturados.

Nossa tarefa é resolver os detalhes e, para cada um deles, determinar a probabilidade de que ele tenha sido produzido na máquina nº 1 ou na máquina nº 2. Também assumiremos que a peça foi produzida na máquina, cuja probabilidade é maior.

Solução

Antes de tudo, analisaremos o algoritmo da solução. Podemos calcular facilmente a probabilidade de cada hipótese, no entanto, já sabemos seus valores em um exemplo anterior:P(B1)=0.6, P(B2)=0.4. Mas como calculamos a probabilidade de apreender um produto padrão individualmente para cada uma das hipóteses disponíveis? Tudo é simples! Não consideraremos a probabilidade; em vez disso, determinaremos a densidade de probabilidade do item a ser removido com seu valor de diâmetro para cada hipótese separadamente.

Para fazer isso, usamos a fórmula conhecida:

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



Onde xi- um valor aleatório (no nosso caso, o tamanho real do diâmetro do produto),

μj- expectativa matemática de variáveis ​​aleatórias jhipótese (no nosso caso, o tamanho médio do diâmetro do produto produzido em jmáquina-ferramenta)

σj- desvio padrão de variáveis ​​aleatórias jhipótese (no nosso caso, o desvio do tamanho médio do diâmetro do produto produzido em jmáquina-ferramenta).

Assim, como fazemos uma substituição inteligente da probabilidade de densidade de probabilidade no numerador da fórmula de Bayes, realizaremos substituições semelhantes no denominador.

Nós adaptamos as fórmulas especificamente para o nosso caso.

A probabilidade de a peça com dimensões de diâmetroxi produzido na máquina nº 1 é determinado pela fórmula:

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



A probabilidade de a peça com dimensões de diâmetro xiproduzido na máquina nº 2:

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



Observe que substituímos a notação pela probabilidade de uma hipótese por P(Bj)no wj. Isso se deve à designação correspondente no curso previamente designado. Além disso, usaremos exatamente essa notação.

Agora estamos 100% completos e prontos para resolver o problema.

Nós olhamos para o código

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


Funções usadas no segundo exemplo
#      ,    
#      : .,  . 
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()


Especifique os parâmetros
#    
#      №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.


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


Gráfico # 1



O que acabamos de fazer? Geramos de maneira pseudo-aleatória 10.000 valores descritos pela distribuição normal, dos quais 6.000 com média de 64 mm e desvio padrão de 4 mm, 4.000 com média de 52 mm e desvio padrão de 2 mm. Além disso, de acordo com o algoritmo acima, para cada parte, foi determinada a probabilidade de sua produção na máquina No. 1 ou No. 2. Depois disso, escolhemos a hipótese sobre a produção do produto em uma ou outra máquina, dependendo da hipótese mais provável. Finalmente, comparamos os resultados do nosso algoritmo com as respostas verdadeiras.

Na saída, obtivemos uma parte das respostas corretas - 0,986. Em geral, isso é muito bom, uma vez que realizamos o treinamento sem usar respostas verdadeiras.

Vamos olhar para o gráfico. O que é recomendado prestar atenção? Veja onde estão localizados os produtos que não são identificados corretamente pelo algoritmo.

Primeiro, vemos erros de algoritmo apenas na área em que produtos fabricados em máquinas diferentes têm o mesmo diâmetro. Parece bastante lógico.

Em segundo lugar, vemos que o algoritmo está errado apenas nos objetos que estão mais distantes do verdadeiro valor médio dos objetos e, ao mesmo tempo, estão bem próximos do centro falso.

Mas o mais interessante é que os problemas começam principalmente depois de cruzar valores aproximadamente iguaisμ+2σ, especificamente no nosso caso no cruzamento μ12σ1e μ2+2σ2

Quem desejar pode "distorcer" os parâmetros e ver como a qualidade do algoritmo muda. Isso será útil para uma melhor assimilação do material.

Bem, passamos ao próximo exemplo

O terceiro exemplo. Caso bidimensional


Desta vez, consideraremos o caso em que temos dados não apenas sobre o diâmetro dos produtos, mas também, por exemplo, sobre o peso de cada produto. Deixe a máquina nº 1 produzir peças com peso de 12 ge um desvio padrão de 1 g, a máquina nº 2 produz produtos com peso de 10 ge um desvio padrão de 0,8 mm.

Para o resto, usaremos os dados do exemplo anterior, anotá-los.

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

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

A situação é a mesma - todos os detalhes são misturados em uma grande pilha e precisamos separá-los e determinar para cada parte a probabilidade de sua produção na máquina nº 1 e na máquina nº 2.

Solução

Em geral, a abordagem da solução do exemplo anterior não muda, com exceção da fórmula para determinar a densidade de probabilidade de uma variável aleatória. Para o caso multidimensional, a seguinte fórmula é usada:

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



Onde ΣjÉ uma matriz de covariâncias de variáveis ​​aleatórias jhipótese (no nosso caso, a matriz de covariâncias dos parâmetros dos produtos produzidos em jmáquina-ferramenta)

ΣjÉ o determinante da matriz de covariância de variáveis ​​aleatórias jhipótese

μj- expectativa matemática de variáveis ​​aleatórias jhipótese (no nosso caso, o valor médio dos produtos fabricados em jmáquina-ferramenta)

xi- variável aleatória (no nosso caso, os parâmetros i-ésimo produto)

Desta vez, para entender o código, um leitor respeitado precisará de conhecimento sobre operações com matrizes

Funções usadas no terceiro exemplo
#      ,    
#      : .,  . 
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()


Especifique os parâmetros
#  
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


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


Gráfico No. 2.1 “Distribuição de produtos de acordo com o algoritmo”


Gráfico No. 2.2 “Distribuição verdadeira de produtos”


Por analogia com o exemplo anterior, geramos 10.000 valores de acordo com os parâmetros acimaμ e σ, escreveu várias funções para a operação do nosso algoritmo e o lançou. A diferença fundamental entre esse código e o código do exemplo anterior é que desta vez usamos expressões matriciais para cálculos.

Como resultado, nosso algoritmo mostra a proporção de respostas corretas igual a 0,998, o que é realmente muito bom.

Os problemas, como vemos, são todos iguais - erros nos detalhes que foram produzidos em máquinas diferentes e ao mesmo tempo têm dimensões e pesos semelhantes.

O código é escrito para que o leitor possa substituir qualquer valor de parâmetroμ e σe veja como o algoritmo funcionará: em quais casos a qualidade se deteriorará e em que melhorará. O principal é não esquecer de estudar a programação. Bem, estamos passando para o nosso ponto final na análise do algoritmo EM, na verdade, o próprio algoritmo EM.

Conheça o algoritmo EM


Continuamos nosso exemplo com as peças que chegaram ao armazém. Mas desta vez saberemos apenas que os produtos foram fabricados em duas máquinas diferentes, existem 10.000 delas, cada peça tem diâmetro e tamanho e não sabemos mais nada. Mas a tarefa não mudou - nós, como antes, de toda a pilha grande, misturamos aleatoriamente produtos, precisaremos determinar a qual máquina essa ou aquela parte pertence.

À primeira vista, parece quase irreal, mas na verdade temos em nossas mãos um poderoso conjunto de ferramentas: a fórmula de Bayes e a fórmula de densidade de probabilidade de uma variável aleatória. Vamos levar toda essa bondade.

Solução

O que faremos? Como deveria ser no algoritmo EM, primeiro inicializamos os parâmetros: A

probabilidade da hipótese é extrair a parte produzida na máquina nº 1 -w1 determinamos a probabilidade de uma hipótese extrair a peça produzida na máquina nº 2 - w2. Existem apenas duas hipóteses, o que significa que cada uma delas na primeira etapa será igual a 0,5.

Expectativa matemática de variáveis ​​aleatóriasμdefina da seguinte maneira. Nós misturamos todos os produtos usando a função aleatória, dividimos o conjunto igualmente em duas partes, para cada parte para cada parâmetro (diâmetro, peso) determinamos o valor médio.

Desvio padrãoσpegue o que é chamado do teto - defina-o como unidade em todos os aspectos. Escrevemos no formato da matriz de covariância.

Estamos prontos para fazer a primeira etapa E do algoritmo. Usando os parâmetros inicializados de variáveis ​​aleatórias, determinamos a probabilidade de cada peça ser atribuída à máquina nº 1 ou à máquina nº 2.

Na verdade, dessa maneira, demos o primeiro passo eletrônico.

Agora cabe ao passo M. Tudo é simples aqui. Depois de determinarmos a probabilidade de cada peça ser produzida em uma máquina específica, podemos recalcular a probabilidade de cada hipótese -w1, w2e μe σ.

Faremos 15 dessas iterações, em duas etapas cada,

e veremos o código.
Funções usadas no algoritmo 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


Especifique os parâmetros
#   
#      №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


Vamos calcular
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()


Quadro 3.1 “Distribuição de peças em cada iteração do algoritmo EM”

Todas as iterações
















O que aconteceu no final? E aconteceu o seguinte, sem informações preliminares sobre a probabilidade de hipóteses para a extração de peças produzidas na máquina nº 1 -w1 ou número da máquina 2 w2, expectativa matemática e desvio padrão de variáveis ​​aleatórias - μe σ, alcançamos uma qualidade comparável ao exemplo anterior, na qual todos os parâmetros acima eram conhecidos.

Compare os parâmetros que definimos ao gerar o conjunto de produtos, os chamados parâmetros verdadeiros e os parâmetros que o algoritmo EM determinou para nós:

Parâmetros verdadeiros

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

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

Parâmetros específicos

μ1164.03,μ1213.99,σ113.44,σ121.23

μ2152.06,μ129.54,σ211.65,σ220.77

Neste artigo chega ao fim. Concluindo, acrescentarei que, para fins de análise de dados, não faz sentido reinventar a roda e escrever o algoritmo EM por conta própria; você pode simplesmente usar a função da biblioteca sklearn fornecida.

Vamos ver como a biblioteca gaussianMixture do sklearn funciona no nosso exemplo.

Módulo de importação
from sklearn.mixture import GaussianMixture as GMM


Especifique os parâmetros
#   
#      №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))


Execute o algoritmo 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)


Trabalhou com um estrondo! Sem erros! Como seria de esperar, a biblioteca GaussianMixture do skussian é mais rápida e melhor que o nosso algoritmo. No entanto, o objetivo do artigo, como lembramos, era entender como os algoritmos EM funcionam. Espero que o objetivo tenha sido alcançado :)

Fontes de informação


Literatura

1) Estatística aplicada: classificações e redução de dimensionalidade, S.A. Ayvazyan, V.M. Buchstaber, I.S. Enyukov, L.D. Meshalkin, Moscou, Finanças e Estatística, 1989

Recursos da Internet

1) Artigo 'Um tutorial delicado do algoritmo EM e sua aplicação à estimativa de parâmetros para modelos gaussianos de mistura e Markov oculto', Instituto Internacional de Ciência da Computação, Jeff A. Bilmes, 1998

2) Apresentação 'Modelos de Mistura e Algoritmo EM', Universidade de NY, David Sontag

3) Curso 'Pesquisa por estruturas de dados'

4) Fórmula para probabilidade total e fórmulas de Bayes

5) Função 'sklearn.mixture.GaussianMixture'

6) Distribuição normal, Wikipedia

7) Multidimensional distribuição normal, Wikipedia

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

9) Um artigo sobre o algoritmo EM

10) Outro artigo sobre o algoritmo EM

11) Github github.com/AlexanderPetrenko83/Articles

All Articles