Числовой пример для понимания максимизации ожидания

Я пытаюсь понять алгоритм EM, чтобы иметь возможность его реализовать и использовать. Я провел целый день, читая теорию и документ, где EM используется для отслеживания самолета с использованием информации о местоположении, поступающей с радара. Честно говоря, я не думаю, что полностью понимаю основную идею. Может кто-нибудь указать мне на числовой пример, показывающий несколько итераций (3-4) EM для более простой задачи (например, оценки параметров гауссовского распределения или последовательности синусоидального ряда или подгонки линии).

Даже если кто-то может указать мне на кусок кода (с синтетическими данными), я могу попытаться пройтись по коду.

Это рецепт изучения EM на практическом и (на мой взгляд) очень интуитивном примере с Coin-Toss:

1. Прочитайте этот короткий учебник по EM от Do и Batzoglou. Это схема, где объясняется пример броска монеты:

   

2. У вас могут быть вопросительные знаки в вашей голове, особенно относительно того, откуда берутся вероятности на этапе Ожидания. Пожалуйста, ознакомьтесь с объяснениями на этой странице обмена стека математики .

3. Посмотрите / запустите этот код, который я написал на Python, который имитирует решение проблемы броска монеты, в учебном материале по EM из пункта 1:

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

Похоже, ваш вопрос состоит из двух частей: основной идеи и конкретного примера. Я начну с основной идеи, затем сошлюсь на пример внизу.

EM полезно в Catch-22 ситуаций , когда кажется , что вы должны знать , прежде чем можно вычислить , и вы должны знать , прежде чем можно вычислить .B B $A$$B$$B$$A$

Наиболее распространенный случай, с которым сталкиваются люди, это, вероятно, смешанные распределения. Для нашего примера давайте рассмотрим простую модель гауссовой смеси:

У вас есть два разных одномерных гауссовых распределения с разными средними и единичной дисперсией.

У вас есть куча точек данных, но вы не уверены, какие точки пришли из какого дистрибутива, и вы также не уверены в способах двух дистрибутивов.

А теперь ты застрял

• Если бы вы знали истинное значение, вы могли бы выяснить, какие точки данных были получены с какой гауссианы. Например, если точка данных имела очень высокое значение, она, вероятно, была получена из распределения с более высоким средним значением. Но вы не знаете, каковы средства, поэтому это не сработает.

• Если вы знали, из какого распределения пришла каждая точка, то вы могли бы оценить средние значения двух распределений, используя выборочные средние значения соответствующих точек. Но вы на самом деле не знаете, какие точки назначить какому распределению, так что это тоже не сработает.

Таким образом, ни один из подходов не выглядит так, как будто он работает: вам нужно знать ответ, прежде чем вы сможете найти ответ, и вы застряли.

Что EM позволяет вам делать, так это чередовать эти два поддающихся обработке шага вместо одновременного выполнения всего процесса.

Вам нужно начать с предположения о двух способах (хотя ваше предположение не обязательно должно быть очень точным, вам нужно с чего-то начинать).

Если ваше предположение о средстве было точным, то у вас было бы достаточно информации, чтобы выполнить шаг, описанный в моем первом пункте выше, и вы могли бы (вероятностно) назначить каждую точку данных одному из двух гауссиан. Хотя мы знаем, что наше предположение неверно, давайте все равно попробуем. И затем, учитывая назначенные распределения каждой точки, вы можете получить новые оценки для средних с использованием второй маркированной точки. Оказывается, что каждый раз, выполняя эти два шага, вы улучшаете нижнюю границу вероятности модели.

Это уже довольно круто: даже несмотря на то, что два предложения в пунктах выше не выглядят так, как будто они будут работать индивидуально, вы все равно можете использовать их вместе, чтобы улучшить модель. Настоящая магия ЕСТ является то, что после достаточно итераций, нижняя граница будет настолько высока , что не будет никакого пространства между ним и локальным максимумом. В результате и вы локально оптимизировали вероятность.

Таким образом, вы не просто улучшили модель, вы нашли лучшую возможную модель, которую можно найти с помощью дополнительных обновлений.

На этой странице из Википедии показан немного более сложный пример (двумерные гауссианы и неизвестная ковариация), но основная идея та же. Он также включает в себя хорошо прокомментированный Rкод для реализации примера.

В коде шаг «Ожидание» (электронный шаг) соответствует моему первому пункту: выяснить, какой гауссов отвечает за каждую точку данных, учитывая текущие параметры для каждого гауссиана. Шаг «Максимизация» (M-шаг) обновляет средние значения и ковариации, учитывая эти назначения, как в моем втором пункте пули.

Как вы можете видеть из анимации, эти обновления быстро позволяют алгоритму перейти от набора ужасных оценок к набору очень хороших: на самом деле, похоже, существуют два облака точек с центром в двух гауссовых распределениях, которые находит EM.

Вот пример максимизации ожидания (EM), используемой для оценки среднего значения и стандартного отклонения. Код написан на Python, но за ним должно быть легко следовать, даже если вы не знакомы с языком.

Мотивация для EM

Красные и синие точки, показанные ниже, взяты из двух разных нормальных распределений, каждое из которых имеет определенное среднее значение и стандартное отклонение:

Чтобы вычислить разумные аппроксимации «истинного» среднего значения и параметров стандартного отклонения для красного распределения, мы могли бы очень легко посмотреть на красные точки и записать положение каждой из них, а затем использовать знакомые формулы (и аналогично для синей группы) ,

Теперь рассмотрим случай, когда мы знаем, что есть две группы точек, но мы не можем видеть, какая точка принадлежит какой группе. Другими словами, цвета скрыты:

Совсем не очевидно, как разделить точки на две группы. Теперь мы не можем просто посмотреть на позиции и вычислить оценки для параметров красного или синего распределения.

Здесь EM может быть использован для решения проблемы.

Использование EM для оценки параметров

Вот код, используемый для генерации точек, показанных выше. Вы можете увидеть фактические средние и стандартные отклонения нормальных распределений, из которых были нарисованы точки. Переменные redи blueсодержат позиции каждой точки в красной и синей группах соответственно:

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

Если бы мы могли видеть цвет каждой точки, мы попытались бы восстановить средние значения и стандартные отклонения, используя библиотечные функции:

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

Но так как цвета скрыты от нас, мы начнем процесс EM ...

Сначала мы просто угадываем значения для параметров каждой группы ( шаг 1 ). Эти догадки не должны быть хорошими:

# 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

Довольно плохие догадки - средства выглядят так, будто они далеки от любой «середины» группы точек.

Чтобы продолжить с EM и улучшить эти догадки, мы вычисляем вероятность того, что каждая точка данных (независимо от ее секретного цвета) появится под этими догадками для среднего и стандартного отклонения ( шаг 2 ).

Переменная both_coloursсодержит каждую точку данных. Функция stats.normвычисляет вероятность точки при нормальном распределении с заданными параметрами:

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)

Это говорит нам, например, что с нашими текущими догадками точка данных в 1.761, скорее всего, будет красной (0.189), чем синей (0.00003).

Мы можем превратить эти два значения правдоподобия в веса ( шаг 3 ), чтобы они суммировали до 1 следующим образом:

likelihood_total = likelihood_of_red + likelihood_of_blue

red_weight = likelihood_of_red / likelihood_total
blue_weight = likelihood_of_blue / likelihood_total

С нашими текущими оценками и нашими недавно вычисленными весами мы можем теперь вычислить новые, возможно, лучшие оценки параметров ( шаг 4 ). Нам нужна функция для среднего значения и функция для стандартного отклонения:

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)

Они выглядят очень похоже на обычные функции для среднего значения и стандартного отклонения данных. Разница заключается в использовании weightпараметра, который присваивает вес каждой точке данных.

Это взвешивание является ключом к EM. Чем больше вес цвета в точке данных, тем больше точка данных влияет на следующие оценки параметров этого цвета. В конечном итоге это приводит к вытягиванию каждого параметра в правильном направлении.

Новые догадки вычисляются с помощью этих функций:

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

Процесс ЭМ затем повторяется с этими новыми догадками, начиная с шага 2 и далее. Мы можем повторять шаги для заданного количества итераций (скажем, 20) или пока мы не увидим, что параметры сходятся.

После пяти итераций мы видим, что наши начальные плохие догадки начинают поправляться:

После 20 итераций процесс EM более или менее сходится:

Для сравнения, вот результаты EM-процесса по сравнению со значениями, вычисленными, когда информация о цвете не скрыта:

          | 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

Примечание: этот ответ был адаптирован из моего ответа о переполнении стека здесь .

После ответа Жубарба я реализовал пример EM и Doz Batzoglou «подбрасывание монет» в GNU R. Обратите внимание, что я использую mleфункцию stats4пакета - это помогло мне более четко понять, как связаны EM и MLE.

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

Все вышеперечисленное выглядит как замечательные ресурсы, но я должен сослаться на этот замечательный пример. Это представляет очень простое объяснение для нахождения параметров для двух линий набора точек. Учебник Яир Вайсс в то время как в Массачусетском технологическом институте.

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

Ответ, данный Жубарбом, великолепен, но, к сожалению, он написан на Python. Ниже приведена Java-реализация EM-алгоритма, выполненного для той же задачи (изложена в статье Do and Batzoglou, 2008). Я добавил несколько printf в стандартный вывод, чтобы увидеть, как сходятся параметры.

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

Код Java следует ниже:

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, jandrader@yahoo.com

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

Ну, я бы посоветовал вам прочитать книгу Марии Риццо о R. Одна из глав содержит использование алгоритма EM с числовым примером. Я помню, как просматривал код для лучшего понимания.

Кроме того, попробуйте просмотреть его с точки зрения кластеризации в начале. Разработайте вручную задачу кластеризации, где 10 наблюдений взяты из двух разных нормальных плотностей. Это должно помочь. Воспользуйтесь помощью R :)

$\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