Numeryczny przykład zrozumienia maksymalizacji oczekiwań

Staram się dobrze zrozumieć algorytm EM, aby móc go wdrożyć i używać. Spędziłem cały dzień czytając teorię i artykuł, w którym EM służy do śledzenia samolotu z wykorzystaniem informacji o położeniu pochodzących z radaru. Szczerze mówiąc, nie sądzę, że w pełni rozumiem leżącą u podstaw ideę. Czy ktoś może wskazać mi numeryczny przykład pokazujący kilka iteracji EM (3–4) dla prostszego problemu (np. Oszacowanie parametrów rozkładu Gaussa lub sekwencji szeregu sinusoidalnego lub dopasowanie linii).

Nawet jeśli ktoś może wskazać mi fragment kodu (z danymi syntetycznymi), mogę spróbować przejść przez kod.

To jest przepis na naukę EM z praktycznym i (moim zdaniem) bardzo intuicyjnym przykładem „Rzut monetą”:

1. Przeczytaj ten krótki samouczek EM autorstwa Do i Batzoglou. Oto schemat, w którym wyjaśniono przykład rzucania monetą:

   

2. Możesz mieć w głowie znaki zapytania, szczególnie jeśli chodzi o to, skąd pochodzą prawdopodobieństwa na etapie Oczekiwania. Proszę spojrzeć na wyjaśnienia na tej stronie wymiany stosów matematycznych .

3. Spójrz na / uruchom ten kod, który napisałem w Pythonie, który symuluje rozwiązanie problemu rzucania monetą w dokumencie instruktażowym EM w punkcie 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()
   

Wygląda na to, że twoje pytanie składa się z dwóch części: idei i konkretnego przykładu. Zacznę od idei leżącej u podstaw, a następnie link do przykładu na dole.

$A$B $B$$B$$A$

Najczęstszym przypadkiem, z którym ludzie mają do czynienia, jest prawdopodobnie rozkład mieszanin. W naszym przykładzie spójrzmy na prosty model mieszanki Gaussa:

Masz dwie różne jednowymiarowe rozkłady Gaussa z różnymi średnimi i wariancjami jednostek.

Masz kilka punktów danych, ale nie masz pewności, które punkty pochodzą z której dystrybucji, a także nie jesteś pewien co do średnich dwóch dystrybucji.

A teraz utknąłeś:

• Jeśli znasz prawdziwe środki, możesz dowiedzieć się, które punkty danych pochodzą z którego Gaussa. Na przykład, jeśli punkt danych miał bardzo wysoką wartość, prawdopodobnie pochodzi on z rozkładu o wyższej średniej. Ale nie wiesz, jakie są środki, więc to nie zadziała.

• Gdybyście wiedzieli, z którego rozkładu pochodzi każdy punkt, moglibyście oszacować te dwa rozkłady za pomocą średnich próbek odpowiednich punktów. Ale tak naprawdę nie wiesz, które punkty przypisać do której dystrybucji, więc to też nie zadziała.

Żadne podejście nie wydaje się więc skuteczne: musisz znać odpowiedź, zanim znajdziesz odpowiedź, i utkniesz.

EM pozwala ci na zmianę na przemian między tymi dwoma wykonalnymi krokami, zamiast zajmować się całym procesem naraz.

Musisz zacząć od zgadywania na temat dwóch środków (chociaż twoje przypuszczenie niekoniecznie musi być bardzo dokładne, musisz gdzieś zacząć).

Jeśli twoje przypuszczenia na temat średnich były dokładne, miałbyś wystarczającą ilość informacji, aby wykonać krok w moim pierwszym punkcie powyżej, i możesz (probabilistycznie) przypisać każdy punkt danych do jednego z dwóch Gaussów. Chociaż wiemy, że nasze domysły są błędne, spróbujmy mimo to. A następnie, biorąc pod uwagę przypisane rozkłady każdego punktu, możesz uzyskać nowe oszacowania średnich za pomocą drugiego punktu. Okazuje się, że za każdym razem, gdy wykonujesz pętlę przez te dwa kroki, poprawiasz dolną granicę prawdopodobieństwa modelu.

To już całkiem fajne: chociaż dwie sugestie w punktach powyżej nie wyglądały tak, jakby działały indywidualnie, nadal możesz je wykorzystać razem, aby ulepszyć model. Prawdziwa magia EM jest to, że po tyle powtórzeń, dolna granica będzie tak wysoka, że nie będzie żadnych przestrzeń między nią a lokalnym maksimum. W rezultacie lokalnie zoptymalizowałeś prawdopodobieństwo.

Więc nie tylko ulepszyłeś model, ale znalazłeś najlepszy możliwy model, jaki można znaleźć dzięki aktualizacjom przyrostowym.

Ta strona z Wikipedii pokazuje nieco bardziej skomplikowany przykład (dwuwymiarowe Gaussianie i nieznana kowariancja), ale podstawowa idea jest taka sama. Zawiera również dobrze skomentowany Rkod do implementacji przykładu.

W kodzie krok „Oczekiwanie” (krok E) odpowiada mojemu pierwszemu punktowi: ustalenie, który Gaussian bierze odpowiedzialność za każdy punkt danych, biorąc pod uwagę bieżące parametry dla każdego Gaussa. Krok „Maksymalizacja” (krok M) aktualizuje środki i kowariancje, biorąc pod uwagę te zadania, jak w moim drugim punkcie.

Jak widać na animacji, te aktualizacje szybko pozwalają algorytmowi przejść od zestawu okropnych oszacowań do zestawu bardzo dobrych: naprawdę wydaje się, że istnieją dwie chmury punktów wyśrodkowane na dwóch rozkładach Gaussa, które znajdzie EM.

Oto przykład metody Expectation Maximization (EM) zastosowanej do oszacowania średniej i odchylenia standardowego. Kod jest w języku Python, ale powinien być łatwy do naśladowania, nawet jeśli nie znasz języka.

Motywacja do EM

Przedstawione poniżej czerwone i niebieskie punkty pochodzą z dwóch różnych rozkładów normalnych, z których każdy ma określoną średnią i standardowe odchylenie:

Aby obliczyć rozsądne aproksymacje „prawdziwych” średnich i standardowych parametrów odchylenia dla rozkładu czerwonego, możemy bardzo łatwo spojrzeć na czerwone punkty i zapisać położenie każdego z nich, a następnie użyć znanych wzorów (i podobnie dla grupy niebieskiej) .

Rozważmy teraz przypadek, w którym wiemy, że istnieją dwie grupy punktów, ale nie widzimy, który punkt należy do której grupy. Innymi słowy, kolory są ukryte:

Nie jest wcale oczywiste, jak podzielić punkty na dwie grupy. Nie jesteśmy teraz w stanie spojrzeć na pozycje i obliczyć oszacowań parametrów rozkładu czerwonego lub niebieskiego.

Tutaj EM można wykorzystać do rozwiązania problemu.

Wykorzystanie EM do oszacowania parametrów

Oto kod używany do generowania punktów pokazanych powyżej. Możesz zobaczyć rzeczywiste średnie i standardowe odchylenia rozkładów normalnych, z których zostały narysowane punkty. Zmienne redi blueutrzymują pozycje każdego punktu odpowiednio w grupie czerwonej i niebieskiej:

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

Gdybyśmy mogli zobaczyć kolor każdego punktu, spróbowalibyśmy odzyskać średnie i standardowe odchylenia za pomocą funkcji bibliotecznych:

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

Ale ponieważ kolory są przed nami ukryte, rozpoczniemy proces EM ...

Po pierwsze, domyślamy się tylko wartości parametrów każdej grupy ( krok 1 ). Te domysły nie muszą być dobre:

# 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

Całkiem złe domysły - środki wyglądają, jakby były daleko od jakiegokolwiek „środka” grupy punktów.

Aby kontynuować EM i poprawić te przypuszczenia, obliczamy prawdopodobieństwo, że każdy punkt danych (niezależnie od jego tajnego koloru) pojawi się pod tymi przypuszczeniami dla średniej i odchylenia standardowego ( krok 2 ).

Zmienna both_coloursprzechowuje każdy punkt danych. Funkcja stats.normoblicza prawdopodobieństwo punktu o rozkładzie normalnym na podstawie podanych parametrów:

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)

To mówi nam na przykład, że przy naszych bieżących przypuszczeniach punkt danych na 1,761 jest znacznie bardziej czerwony (0,189) niż niebieski (0,00003).

Możemy zamienić te dwie wartości prawdopodobieństwa na wagi ( krok 3 ), aby sumowały się do 1 w następujący sposób:

likelihood_total = likelihood_of_red + likelihood_of_blue

red_weight = likelihood_of_red / likelihood_total
blue_weight = likelihood_of_blue / likelihood_total

Dzięki naszym bieżącym oszacowaniom i naszym nowo obliczonym wagom możemy teraz obliczyć nowe, prawdopodobnie lepsze oszacowania parametrów ( krok 4 ). Potrzebujemy funkcji dla średniej i funkcji dla odchylenia standardowego:

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)

Wyglądają one bardzo podobnie do zwykłych funkcji do średniej i odchylenia standardowego danych. Różnica polega na zastosowaniu weightparametru, który przypisuje wagę do każdego punktu danych.

Ta waga jest kluczem do EM. Im większa waga koloru w punkcie danych, tym bardziej punkt danych wpływa na następne oszacowania parametrów tego koloru. Ostatecznie powoduje to pociągnięcie każdego parametru we właściwym kierunku.

Nowe domysły są obliczane za pomocą tych funkcji:

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

Proces EM jest następnie powtarzany z tymi nowymi domysłami, począwszy od kroku 2. Możemy powtórzyć kroki dla danej liczby iteracji (powiedzmy 20) lub dopóki nie zobaczymy, że parametry się zbiegają.

Po pięciu iteracjach widzimy, że nasze początkowe błędne domysły zaczynają się poprawiać:

Po 20 iteracjach proces EM zbliżył się mniej więcej:

Dla porównania, oto wyniki procesu EM w porównaniu z wartościami obliczonymi, gdy informacja o kolorze nie jest ukryta:

          | 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

Uwaga: ta odpowiedź została zaadaptowana z mojej odpowiedzi na temat Przepełnienia stosu tutaj .

Po odpowiedzi Zhubarba zaimplementowałem przykład EM „Do podrzucania monet” Do i Batzoglou w GNU R. Zauważ, że korzystam z mlefunkcji stats4pakietu - pomogło mi to lepiej zrozumieć związek EM i 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))

Wszystkie powyższe wyglądają jak świetne zasoby, ale muszę odnieść się do tego wspaniałego przykładu. Przedstawia bardzo proste wyjaśnienie dotyczące znalezienia parametrów dla dwóch linii zestawu punktów. Samouczek przygotował Yair Weiss podczas MIT.

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

Odpowiedź udzielona przez Zhubarb jest świetna, ale niestety jest w Pythonie. Poniżej znajduje się implementacja Java algorytmu EM wykonana na tym samym problemie (postawiona w artykule Do i Batzoglou, 2008). Dodałem kilka printf do standardowego wyjścia, aby zobaczyć, jak parametry się zbiegają.

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

Kod Java poniżej:

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

Sugerowałbym, żebyś przejrzał książkę Marii L. Rizzo na temat R. Jeden z rozdziałów zawiera użycie algorytmu EM z przykładem numerycznym. Pamiętam, jak przechodziłem przez kod dla lepszego zrozumienia.

Spróbuj też wyświetlić go na początku z punktu widzenia grupowania. Opracuj ręcznie, problem grupowania, w którym 10 obserwacji jest pobieranych z dwóch różnych normalnych gęstości. To powinno pomóc. Weź pomoc od 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