Exemplo numérico para entender a Expectativa-Maximização

Estou tentando entender bem o algoritmo EM, para poder implementá-lo e usá-lo. Passei um dia inteiro lendo a teoria e um artigo em que o EM é usado para rastrear uma aeronave usando as informações de posição provenientes de um radar. Honestamente, acho que não entendo completamente a idéia subjacente. Alguém pode me apontar para um exemplo numérico que mostra algumas iterações (3-4) do EM para um problema mais simples (como estimar os parâmetros de uma distribuição gaussiana ou uma sequência de uma série sinusoidal ou ajustar uma linha).

Mesmo que alguém possa me indicar um pedaço de código (com dados sintéticos), eu posso tentar percorrer o código.

Esta é uma receita para aprender EM com um exemplo prático e (na minha opinião) muito intuitivo 'Coin-Toss':

1. Leia este pequeno tutorial sobre EM de Do e Batzoglou. Este é o esquema no qual o exemplo do sorteio é explicado:

   

2. Você pode ter pontos de interrogação em sua mente, principalmente em relação à origem das probabilidades na etapa Expectativa. Veja as explicações nesta página de troca de pilhas matemáticas .

3. Veja / execute este código que escrevi em Python que simula a solução para o problema do sorteio no artigo do tutorial 1 do EM:

   import numpy as np
   import math
   import matplotlib.pyplot as plt
   
   ## E-M Coin Toss Example as given in the EM tutorial paper by Do and Batzoglou* ##
   
   def get_binomial_log_likelihood(obs,probs):
       """ Return the (log)likelihood of obs, given the probs"""
       # Binomial Distribution Log PDF
       # ln (pdf)      = Binomial Coeff * product of probabilities
       # ln[f(x|n, p)] =   comb(N,k)    * num_heads*ln(pH) + (N-num_heads) * ln(1-pH)
   
       N = sum(obs);#number of trials  
       k = obs[0] # number of heads
       binomial_coeff = math.factorial(N) / (math.factorial(N-k) * math.factorial(k))
       prod_probs = obs[0]*math.log(probs[0]) + obs[1]*math.log(1-probs[0])
       log_lik = binomial_coeff + prod_probs
   
       return log_lik
   
   # 1st:  Coin B, {HTTTHHTHTH}, 5H,5T
   # 2nd:  Coin A, {HHHHTHHHHH}, 9H,1T
   # 3rd:  Coin A, {HTHHHHHTHH}, 8H,2T
   # 4th:  Coin B, {HTHTTTHHTT}, 4H,6T
   # 5th:  Coin A, {THHHTHHHTH}, 7H,3T
   # so, from MLE: pA(heads) = 0.80 and pB(heads)=0.45
   
   # represent the experiments
   head_counts = np.array([5,9,8,4,7])
   tail_counts = 10-head_counts
   experiments = zip(head_counts,tail_counts)
   
   # initialise the pA(heads) and pB(heads)
   pA_heads = np.zeros(100); pA_heads[0] = 0.60
   pB_heads = np.zeros(100); pB_heads[0] = 0.50
   
   # E-M begins!
   delta = 0.001  
   j = 0 # iteration counter
   improvement = float('inf')
   while (improvement>delta):
       expectation_A = np.zeros((len(experiments),2), dtype=float) 
       expectation_B = np.zeros((len(experiments),2), dtype=float)
       for i in range(0,len(experiments)):
           e = experiments[i] # i'th experiment
             # loglikelihood of e given coin A:
           ll_A = get_binomial_log_likelihood(e,np.array([pA_heads[j],1-pA_heads[j]])) 
             # loglikelihood of e given coin B
           ll_B = get_binomial_log_likelihood(e,np.array([pB_heads[j],1-pB_heads[j]])) 
   
             # corresponding weight of A proportional to likelihood of A 
           weightA = math.exp(ll_A) / ( math.exp(ll_A) + math.exp(ll_B) ) 
   
             # corresponding weight of B proportional to likelihood of B
           weightB = math.exp(ll_B) / ( math.exp(ll_A) + math.exp(ll_B) ) 
   
           expectation_A[i] = np.dot(weightA, e) 
           expectation_B[i] = np.dot(weightB, e)
   
       pA_heads[j+1] = sum(expectation_A)[0] / sum(sum(expectation_A)); 
       pB_heads[j+1] = sum(expectation_B)[0] / sum(sum(expectation_B)); 
   
       improvement = ( max( abs(np.array([pA_heads[j+1],pB_heads[j+1]]) - 
                       np.array([pA_heads[j],pB_heads[j]]) )) )
       j = j+1
   
   plt.figure();
   plt.plot(range(0,j),pA_heads[0:j], 'r--')
   plt.plot(range(0,j),pB_heads[0:j])
   plt.show()
   

Parece que sua pergunta tem duas partes: a ideia subjacente e um exemplo concreto. Começarei com a ideia subjacente e, em seguida, vinculo a um exemplo na parte inferior.

EM é útil em Catch-22 situações em que parece que você precisa saber antes de poder calcular e você precisa saber antes de poder calcular .B B $A$$B$$B$$A$

O caso mais comum com o qual as pessoas lidam é provavelmente a distribuição de misturas. Para o nosso exemplo, vejamos um modelo simples de mistura gaussiana:

Você tem duas distribuições Gaussianas univariadas diferentes, com diferentes médias e variação de unidades.

Você tem vários pontos de dados, mas não tem certeza de quais pontos vieram de qual distribuição e também não tem certeza sobre os meios das duas distribuições.

E agora você está preso:

• Se você soubesse o verdadeiro meio, poderia descobrir quais pontos de dados vieram de quais gaussianos. Por exemplo, se um ponto de dados tinha um valor muito alto, provavelmente vinha da distribuição com a média mais alta. Mas você não sabe quais são os meios, então isso não vai funcionar.

• Se você soubesse de qual distribuição cada ponto veio, poderia estimar as médias das duas distribuições usando as médias amostrais dos pontos relevantes. Mas você realmente não sabe quais pontos atribuir a qual distribuição, portanto, isso também não funcionará.

Portanto, nenhuma das abordagens parece funcionar: você precisaria saber a resposta antes de poder encontrá-la e ficar preso.

O que o EM permite fazer é alternar entre essas duas etapas tratáveis, em vez de enfrentar todo o processo de uma só vez.

Você precisará começar com um palpite sobre os dois meios (embora seu palpite não precise necessariamente ser muito preciso, você precisa começar em algum lugar).

Se o seu palpite sobre os meios fosse preciso, você teria informações suficientes para executar a etapa no meu primeiro ponto acima e poderia (probabilisticamente) atribuir cada ponto de dados a um dos dois gaussianos. Mesmo sabendo que nosso palpite está errado, vamos tentar assim mesmo. E então, dadas as distribuições atribuídas a cada ponto, você pode obter novas estimativas para as médias usando o segundo ponto. Acontece que, cada vez que você percorre essas duas etapas, você está melhorando um limite inferior na probabilidade do modelo.

Isso já é bem legal: mesmo que as duas sugestões nos pontos acima não pareçam funcionar individualmente, você ainda pode usá-las juntas para melhorar o modelo. A verdadeira mágica do EM é que, após iterações suficientes, o limite inferior será tão alto que não haverá espaço entre ele e o máximo local. Como resultado, você otimizou localmente a probabilidade.

Portanto, você não apenas melhorou o modelo, mas encontrou o melhor modelo possível com atualizações incrementais.

Esta página da Wikipedia mostra um exemplo um pouco mais complicado (gaussianos bidimensionais e covariância desconhecida), mas a idéia básica é a mesma. Também inclui Rcódigo bem comentado para implementar o exemplo.

No código, a etapa "Expectativa" (etapa E) corresponde ao meu primeiro ponto: descobrir qual gaussiano é responsável por cada ponto de dados, dados os parâmetros atuais de cada gaussiano. A etapa "Maximização" (etapa M) atualiza os meios e covariâncias, dadas essas atribuições, como no meu segundo item.

Como você pode ver na animação, essas atualizações permitem rapidamente que o algoritmo passe de um conjunto de estimativas terríveis para um conjunto de estimativas muito boas: realmente parece haver duas nuvens de pontos centradas nas duas distribuições gaussianas encontradas pela EM.

Aqui está um exemplo de Maximização de Expectativa (EM) usada para estimar a média e o desvio padrão. O código está em Python, mas deve ser fácil de seguir, mesmo que você não esteja familiarizado com o idioma.

A motivação para EM

Os pontos vermelho e azul mostrados abaixo são extraídos de duas distribuições normais diferentes, cada uma com uma média e desvio padrão específicos:

Para calcular aproximações razoáveis ​​dos parâmetros "verdadeiro" de média e desvio padrão para a distribuição vermelha, poderíamos facilmente olhar para os pontos vermelhos e registrar a posição de cada um e, em seguida, usar as fórmulas familiares (e da mesma forma para o grupo azul) .

Agora considere o caso em que sabemos que existem dois grupos de pontos, mas não podemos ver qual ponto pertence a qual grupo. Em outras palavras, as cores estão ocultas:

Não é de todo óbvio como dividir os pontos em dois grupos. Agora somos incapazes de olhar apenas as posições e calcular estimativas para os parâmetros da distribuição vermelha ou azul.

É aqui que o EM pode ser usado para resolver o problema.

Usando EM para estimar parâmetros

Aqui está o código usado para gerar os pontos mostrados acima. Você pode ver as médias reais e os desvios padrão das distribuições normais das quais os pontos foram extraídos. As variáveis rede bluemantêm as posições de cada ponto nos grupos vermelho e azul respectivamente:

import numpy as np
from scipy import stats

np.random.seed(110) # for reproducible random results

# set parameters
red_mean = 3
red_std = 0.8

blue_mean = 7
blue_std = 2

# draw 20 samples from normal distributions with red/blue parameters
red = np.random.normal(red_mean, red_std, size=20)
blue = np.random.normal(blue_mean, blue_std, size=20)

both_colours = np.sort(np.concatenate((red, blue)))

Se pudéssemos ver a cor de cada ponto, tentaríamos recuperar médias e desvios padrão usando as funções da biblioteca:

>>> np.mean(red)
2.802
>>> np.std(red)
0.871
>>> np.mean(blue)
6.932
>>> np.std(blue)
2.195

Mas como as cores estão ocultas, iniciaremos o processo EM ...

Primeiro, adivinhamos os valores para os parâmetros de cada grupo ( etapa 1 ). Essas suposições não precisam ser boas:

# estimates for the mean
red_mean_guess = 1.1
blue_mean_guess = 9

# estimates for the standard deviation
red_std_guess = 2
blue_std_guess = 1.7

Suposições muito ruins - os meios parecem estar muito longe de qualquer "meio" de um grupo de pontos.

Para continuar com o EM e melhorar essas suposições, calculamos a probabilidade de cada ponto de dados (independentemente de sua cor secreta) aparecer sob essas suposições para a média e o desvio padrão ( etapa 2 ).

A variável both_colourscontém cada ponto de dados. A função stats.normcalcula a probabilidade do ponto em uma distribuição normal com os parâmetros fornecidos:

likelihood_of_red = stats.norm(red_mean_guess, red_std_guess).pdf(both_colours)
likelihood_of_blue = stats.norm(blue_mean_guess, blue_std_guess).pdf(both_colours)

Isso nos diz, por exemplo, que, com nossas suposições atuais, o ponto de dados em 1.761 tem muito mais probabilidade de ser vermelho (0,189) do que azul (0,00003).

Podemos transformar esses dois valores de probabilidade em pesos ( etapa 3 ) para que eles somarem 1 como a seguir:

likelihood_total = likelihood_of_red + likelihood_of_blue

red_weight = likelihood_of_red / likelihood_total
blue_weight = likelihood_of_blue / likelihood_total

Com nossas estimativas atuais e nossos pesos recém-calculados, agora podemos calcular novas estimativas, provavelmente melhores, para os parâmetros ( etapa 4 ). Precisamos de uma função para a média e uma função para o desvio padrão:

def estimate_mean(data, weight):
    return np.sum(data * weight) / np.sum(weight)

def estimate_std(data, weight, mean):
    variance = np.sum(weight * (data - mean)**2) / np.sum(weight)
    return np.sqrt(variance)

Eles parecem muito semelhantes às funções usuais ao desvio médio e padrão dos dados. A diferença é o uso de um weightparâmetro que atribui um peso a cada ponto de dados.

Essa ponderação é a chave para o EM. Quanto maior o peso de uma cor em um ponto de dados, mais ele influencia as próximas estimativas para os parâmetros dessa cor. Por fim, isso tem o efeito de puxar cada parâmetro na direção certa.

As novas suposições são computadas com estas funções:

# new estimates for standard deviation
blue_std_guess = estimate_std(both_colours, blue_weight, blue_mean_guess)
red_std_guess = estimate_std(both_colours, red_weight, red_mean_guess)

# new estimates for mean
red_mean_guess = estimate_mean(both_colours, red_weight)
blue_mean_guess = estimate_mean(both_colours, blue_weight)

O processo EM é então repetido com essas novas suposições da etapa 2 em diante. Podemos repetir as etapas para um determinado número de iterações (digamos 20), ou até vermos os parâmetros convergirem.

Após cinco iterações, vemos nossas primeiras suposições ruins começarem a melhorar:

Após 20 iterações, o processo EM convergiu mais ou menos:

Para comparação, eis os resultados do processo EM comparados com os valores calculados onde as informações de cores não estão ocultas:

          | EM guess | Actual 
----------+----------+--------
Red mean  |    2.910 |   2.802
Red std   |    0.854 |   0.871
Blue mean |    6.838 |   6.932
Blue std  |    2.227 |   2.195

Nota: esta resposta foi adaptada da minha resposta no Stack Overflow aqui .

Seguindo a resposta de Zhubarb, implementei o exemplo de EM de "lançamento de moeda" de Do e Batzoglou no GNU R. Observe que eu uso a mlefunção do stats4pacote - isso me ajudou a entender mais claramente como EM e MLE estão relacionados.

require("stats4");

## sample data from Do and Batzoglou
ds<-data.frame(heads=c(5,9,8,4,7),n=c(10,10,10,10,10),
    coin=c("B","A","A","B","A"),weight_A=1:5*0)

## "baby likelihood" for a single observation
llf <- function(heads, n, theta) {
  comb <- function(n, x) { #nCr function
    return(factorial(n) / (factorial(x) * factorial(n-x)))
  }
  if (theta<0 || theta >1) { # probabilities should be in [0,1]
    return(-Inf);
  }
  z<-comb(n,heads)* theta^heads * (1-theta)^(n-heads);
  return (log(z))
}

## the "E-M" likelihood function
em <- function(theta_A,theta_B) {
  # expectation step: given current parameters, what is the likelihood
  # an observation is the result of tossing coin A (vs coin B)?
  ds$weight_A <<- by(ds, 1:nrow(ds), function(row) {
    llf_A <- llf(row$heads,row$n, theta_A);
    llf_B <- llf(row$heads,row$n, theta_B);

    return(exp(llf_A)/(exp(llf_A)+exp(llf_B)));
  })

  # maximisation step: given params and weights, calculate likelihood of the sample
  return(- sum(by(ds, 1:nrow(ds), function(row) {
    llf_A <- llf(row$heads,row$n, theta_A);
    llf_B <- llf(row$heads,row$n, theta_B);

    return(row$weight_A*llf_A + (1-row$weight_A)*llf_B);
  })))
}

est<-mle(em,start = list(theta_A=0.6,theta_B=0.5), nobs=NROW(ds))

Todas as opções acima parecem ótimos recursos, mas devo vincular a este ótimo exemplo. Apresenta uma explicação muito simples para encontrar os parâmetros para duas linhas de um conjunto de pontos. O tutorial é de Yair Weiss no MIT.

http://www.cs.huji.ac.il/~yweiss/emTutorial.pdf
http://www.cs.huji.ac.il/~yweiss/tutorials.html

A resposta dada por Zhubarb é ótima, mas infelizmente está em Python. Abaixo está uma implementação em Java do algoritmo EM executado no mesmo problema (apresentado no artigo de Do e Batzoglou, 2008). Adicionei alguns printfs à saída padrão para ver como os parâmetros convergem.

thetaA = 0.71301, thetaB = 0.58134
thetaA = 0.74529, thetaB = 0.56926
thetaA = 0.76810, thetaB = 0.54954
thetaA = 0.78316, thetaB = 0.53462
thetaA = 0.79106, thetaB = 0.52628
thetaA = 0.79453, thetaB = 0.52239
thetaA = 0.79593, thetaB = 0.52073
thetaA = 0.79647, thetaB = 0.52005
thetaA = 0.79667, thetaB = 0.51977
thetaA = 0.79674, thetaB = 0.51966
thetaA = 0.79677, thetaB = 0.51961
thetaA = 0.79678, thetaB = 0.51960
thetaA = 0.79679, thetaB = 0.51959
Final result:
thetaA = 0.79678, thetaB = 0.51960

O código Java segue abaixo:

import java.util.*;

/*****************************************************************************
This class encapsulates the parameters of the problem. For this problem posed
in the article by (Do and Batzoglou, 2008), the parameters are thetaA and
thetaB, the probability of a coin coming up heads for the two coins A and B.
*****************************************************************************/
class Parameters
{
    double _thetaA = 0.0; // Probability of heads for coin A.
    double _thetaB = 0.0; // Probability of heads for coin B.

    double _delta = 0.00001;

    public Parameters(double thetaA, double thetaB)
    {
        _thetaA = thetaA;
        _thetaB = thetaB;
    }

    /*************************************************************************
    Returns true if this parameter is close enough to another parameter
    (typically the estimated parameter coming from the maximization step).
    *************************************************************************/
    public boolean converged(Parameters other)
    {
        if (Math.abs(_thetaA - other._thetaA) < _delta &&
            Math.abs(_thetaB - other._thetaB) < _delta)
        {
            return true;
        }

        return false;
    }

    public double getThetaA()
    {
        return _thetaA;
    }

    public double getThetaB()
    {
        return _thetaB;
    }

    public String toString()
    {
        return String.format("thetaA = %.5f, thetaB = %.5f", _thetaA, _thetaB);
    }

}


/*****************************************************************************
This class encapsulates an observation, that is the number of heads
and tails in a trial. The observation can be either (1) one of the
observed observations, or (2) an estimated observation resulting from
the expectation step.
*****************************************************************************/
class Observation
{
    double _numHeads = 0;
    double _numTails = 0;

    public Observation(String s)
    {
        for (int i = 0; i < s.length(); i++)
        {
            char c = s.charAt(i);

            if (c == 'H')
            {
                _numHeads++;
            }
            else if (c == 'T')
            {
                _numTails++;
            }
            else
            {
                throw new RuntimeException("Unknown character: " + c);
            }
        }
    }

    public Observation(double numHeads, double numTails)
    {
        _numHeads = numHeads;
        _numTails = numTails;
    }

    public double getNumHeads()
    {
        return _numHeads;
    }

    public double getNumTails()
    {
        return _numTails;
    }

    public String toString()
    {
        return String.format("heads: %.1f, tails: %.1f", _numHeads, _numTails);
    }

}

/*****************************************************************************
This class runs expectation-maximization for the problem posed by the article
from (Do and Batzoglou, 2008).
*****************************************************************************/
public class EM
{
    // Current estimated parameters.
    private Parameters _parameters;

    // Observations from the trials. These observations are set once.
    private final List<Observation> _observations;

    // Estimated observations per coin. These observations are the output
    // of the expectation step.
    private List<Observation> _expectedObservationsForCoinA;
    private List<Observation> _expectedObservationsForCoinB;

    private static java.io.PrintStream o = System.out;

    /*************************************************************************
    Principal constructor.
    @param observations The observations from the trial.
    @param parameters The initial guessed parameters.
    *************************************************************************/
    public EM(List<Observation> observations, Parameters parameters)
    {
        _observations = observations;
        _parameters = parameters;
    }

    /*************************************************************************
    Run EM until parameters converge.
    *************************************************************************/
    public Parameters run()
    {

        while (true)
        {
            expectation();

            Parameters estimatedParameters = maximization();

            o.printf("%s\n", estimatedParameters);

            if (_parameters.converged(estimatedParameters)) {
                break;
            }

            _parameters = estimatedParameters;
        }

        return _parameters;

    }

    /*************************************************************************
    Given the observations and current estimated parameters, compute new
    estimated completions (distribution over the classes) and observations.
    *************************************************************************/
    private void expectation()
    {

        _expectedObservationsForCoinA = new ArrayList<Observation>();
        _expectedObservationsForCoinB = new ArrayList<Observation>();

        for (Observation observation : _observations)
        {
            int numHeads = (int)observation.getNumHeads();
            int numTails = (int)observation.getNumTails();

            double probabilityOfObservationForCoinA=
                binomialProbability(10, numHeads, _parameters.getThetaA());

            double probabilityOfObservationForCoinB=
                binomialProbability(10, numHeads, _parameters.getThetaB());

            double normalizer = probabilityOfObservationForCoinA +
                                probabilityOfObservationForCoinB;

            // Compute the completions for coin A and B (i.e. the probability
            // distribution of the two classes, summed to 1.0).

            double completionCoinA = probabilityOfObservationForCoinA /
                                     normalizer;
            double completionCoinB = probabilityOfObservationForCoinB /
                                     normalizer;

            // Compute new expected observations for the two coins.

            Observation expectedObservationForCoinA =
                new Observation(numHeads * completionCoinA,
                                numTails * completionCoinA);

            Observation expectedObservationForCoinB =
                new Observation(numHeads * completionCoinB,
                                numTails * completionCoinB);

            _expectedObservationsForCoinA.add(expectedObservationForCoinA);
            _expectedObservationsForCoinB.add(expectedObservationForCoinB);
        }
    }

    /*************************************************************************
    Given new estimated observations, compute new estimated parameters.
    *************************************************************************/
    private Parameters maximization()
    {

        double sumCoinAHeads = 0.0;
        double sumCoinATails = 0.0;
        double sumCoinBHeads = 0.0;
        double sumCoinBTails = 0.0;

        for (Observation observation : _expectedObservationsForCoinA)
        {
            sumCoinAHeads += observation.getNumHeads();
            sumCoinATails += observation.getNumTails();
        }

        for (Observation observation : _expectedObservationsForCoinB)
        {
            sumCoinBHeads += observation.getNumHeads();
            sumCoinBTails += observation.getNumTails();
        }

        return new Parameters(sumCoinAHeads / (sumCoinAHeads + sumCoinATails),
                              sumCoinBHeads / (sumCoinBHeads + sumCoinBTails));

        //o.printf("parameters: %s\n", _parameters);

    }

    /*************************************************************************
    Since the coin-toss experiment posed in this article is a Bernoulli trial,
    use a binomial probability Pr(X=k; n,p) = (n choose k) * p^k * (1-p)^(n-k).
    *************************************************************************/
    private static double binomialProbability(int n, int k, double p)
    {
        double q = 1.0 - p;
        return nChooseK(n, k) * Math.pow(p, k) * Math.pow(q, n-k);
    }

    private static long nChooseK(int n, int k)
    {
        long numerator = 1;

        for (int i = 0; i < k; i++)
        {
            numerator = numerator * n;
            n--;
        }

        long denominator = factorial(k);

        return (long)(numerator / denominator);
    }

    private static long factorial(int n)
    {
        long result = 1;
        for (; n >0; n--)
        {
            result = result * n;
        }

        return result;
    }

    /*************************************************************************
    Entry point into the program.
    *************************************************************************/
    public static void main(String argv[])
    {
        // Create the observations and initial parameter guess
        // from the (Do and Batzoglou, 2008) article.

        List<Observation> observations = new ArrayList<Observation>();
        observations.add(new Observation("HTTTHHTHTH"));
        observations.add(new Observation("HHHHTHHHHH"));
        observations.add(new Observation("HTHHHHHTHH"));
        observations.add(new Observation("HTHTTTHHTT"));
        observations.add(new Observation("THHHTHHHTH"));

        Parameters initialParameters = new Parameters(0.6, 0.5);

        EM em = new EM(observations, initialParameters);

        Parameters finalParameters = em.run();

        o.printf("Final result:\n%s\n", finalParameters);
    }
}

% Implementation of the EM (Expectation-Maximization)algorithm example exposed on:
% Motion Segmentation using EM - a short tutorial, Yair Weiss, %http://www.cs.huji.ac.il/~yweiss/emTutorial.pdf
% Juan Andrade, [email protected]

clear all
clc

%% Setup parameters
m1 = 2;                 % slope line 1
m2 = 6;                 % slope line 2
b1 = 3;                 % vertical crossing line 1
b2 = -2;                % vertical crossing line 2
x = [-1:0.1:5];         % x axis values
sigma1 = 1;             % Standard Deviation of Noise added to line 1
sigma2 = 2;             % Standard Deviation of Noise added to line 2

%% Clean lines
l1 = m1*x+b1;           % line 1
l2 = m2*x+b2;           % line 2

%% Adding noise to lines
p1 = l1 + sigma1*randn(size(l1));
p2 = l2 + sigma2*randn(size(l2));

%% showing ideal and noise values
figure,plot(x,l1,'r'),hold,plot(x,l2,'b'), plot(x,p1,'r.'),plot(x,p2,'b.'),grid

%% initial guess
m11(1) = -1;            % slope line 1
m22(1) = 1;             % slope line 2
b11(1) = 2;             % vertical crossing line 1
b22(1) = 2;             % vertical crossing line 2

%% EM algorithm loop
iterations = 10;        % number of iterations (a stop based on a threshold may used too)

for i=1:iterations

    %% expectation step (equations 2 and 3)
    res1 = m11(i)*x + b11(i) - p1;
    res2 = m22(i)*x + b22(i) - p2;
    % line 1
    w1 = (exp((-res1.^2)./sigma1))./((exp((-res1.^2)./sigma1)) + (exp((-res2.^2)./sigma2)));

    % line 2
    w2 = (exp((-res2.^2)./sigma2))./((exp((-res1.^2)./sigma1)) + (exp((-res2.^2)./sigma2)));

    %% maximization step  (equation 4)
    % line 1
    A(1,1) = sum(w1.*(x.^2));
    A(1,2) = sum(w1.*x);
    A(2,1) = sum(w1.*x);
    A(2,2) = sum(w1);
    bb = [sum(w1.*x.*p1) ; sum(w1.*p1)];
    temp = A\bb;
    m11(i+1) = temp(1);
    b11(i+1) = temp(2);

    % line 2
    A(1,1) = sum(w2.*(x.^2));
    A(1,2) = sum(w2.*x);
    A(2,1) = sum(w2.*x);
    A(2,2) = sum(w2);
    bb = [sum(w2.*x.*p2) ; sum(w2.*p2)];
    temp = A\bb;
    m22(i+1) = temp(1);
    b22(i+1) = temp(2);

    %% plotting evolution of results
    l1temp = m11(i+1)*x+b11(i+1);
    l2temp = m22(i+1)*x+b22(i+1);
    figure,plot(x,l1temp,'r'),hold,plot(x,l2temp,'b'), plot(x,p1,'r.'),plot(x,p2,'b.'),grid
end

Bem, eu sugiro que você leia um livro sobre R de Maria L. Rizzo. Um dos capítulos contém o uso do algoritmo EM com um exemplo numérico. Lembro-me de passar pelo código para entender melhor.

Além disso, tente visualizá-lo do ponto de vista do cluster no começo. Elabore manualmente, um problema de agrupamento em que 10 observações são obtidas de duas densidades normais diferentes. Isso ajuda.Take a ajuda de R :)

Apenas no caso, eu escrevi uma implementação Ruby do exemplo de lançamento de moeda mencionado acima por Do & Batzoglou e produz exatamente os mesmos números que os mesmos parâmetros de entrada e . θ B = $\theta_A = 0.6$$\theta_B = 0.5$

# gem install distribution
require 'distribution'

# error bound
EPS = 10**-6

# number of coin tosses
N = 10

# observations
X = [5, 9, 8, 4, 7]

# randomly initialized thetas
theta_a, theta_b = 0.6, 0.5

p [theta_a, theta_b]

loop do
  expectation = X.map do |h|
    like_a = Distribution::Binomial.pdf(h, N, theta_a)
    like_b = Distribution::Binomial.pdf(h, N, theta_b)

    norm_a = like_a / (like_a + like_b)
    norm_b = like_b / (like_a + like_b)

    [norm_a, norm_b, h]
  end

  maximization = expectation.each_with_object([0.0, 0.0, 0.0, 0.0]) do |(norm_a, norm_b, h), r|
    r[0] += norm_a * h; r[1] += norm_a * (N - h)
    r[2] += norm_b * h; r[3] += norm_b * (N - h)
  end

  theta_a_hat = maximization[0] / (maximization[0] + maximization[1])
  theta_b_hat = maximization[2] / (maximization[2] + maximization[3])

  error_a = (theta_a_hat - theta_a).abs / theta_a
  error_b = (theta_b_hat - theta_b).abs / theta_b

  theta_a, theta_b = theta_a_hat, theta_b_hat

  p [theta_a, theta_b]

  break if error_a < EPS && error_b < EPS
end