EM jest algorytmem maksymalizacji funkcji prawdopodobieństwa, gdy niektóre zmienne w modelu są nieobserwowane (tj. gdy masz zmienne Ukryte).

Można zapytać, jeśli próbujemy zmaksymalizować funkcję, dlaczego nie użyjemy istniejącej maszyny do zmaksymalizowania funkcji. Cóż, jeśli spróbujesz to zmaksymalizować, biorąc pochodne i ustawiając je na zero, okaże się, że w wielu przypadkach warunki pierwszego rzędu nie mają rozwiązania. Jest kurczak z jajkiem problem w tym, aby rozwiązać parametry modelu, musisz znać rozkład nieobserwowanych danych; ale Dystrybucja nieobserwowanych danych jest funkcją parametrów modelu.

E-M próbuje obejść to, iteracyjnie zgadując rozkład nieobserwowanych danych, a następnie szacując parametry modelu, maksymalizując coś, co jest dolną granicą rzeczywistej funkcji prawdopodobieństwa i powtarzając do zbieżności:

Algorytm EM

Start z guess dla wartości parametrów modelu

E-krok: dla każdego punktu danych, który ma brakujące wartości, użyj równania modelu do rozwiązania rozkładu brakujących danych, biorąc pod uwagę bieżące odgadnięcie parametrów modelu i dane obserwowane (zauważ, że rozwiązujesz dla rozkładu dla każdej brakującej wartości, a nie dla wartości oczekiwanej). Teraz, gdy mamy rozkład dla każdej brakującej wartości, możemy obliczyć oczekiwanie funkcji prawdopodobieństwa z szacunek dla nieobserwowanych zmiennych. Jeśli nasze odgadnięcie parametru modelu było poprawne, to oczekiwane prawdopodobieństwo będzie rzeczywistym prawdopodobieństwem naszych obserwowanych danych; jeśli parametry nie były poprawne, będzie to po prostu dolna granica.

M-step: teraz, gdy mamy funkcję oczekiwanego prawdopodobieństwa bez nieobserwowanych zmiennych, zmaksymalizuj funkcję tak, jak w całkowicie obserwowanym przypadku, aby uzyskać nowe oszacowanie parametrów modelu.

Powtarzać aż do konwergencji.

Oto prosty przepis na zrozumienie algorytmu maksymalizacji oczekiwań:

1- przeczytaj ten artykuł em tutorial paper autorstwa Do I Batzoglou.

2- możesz mieć znaki zapytania w głowie, spójrz na wyjaśnienia na tej stronie wymiany stosu matematyki .

3- spójrz na ten kod, który napisałem w Pythonie, który wyjaśnia przykład w artykule z samouczka EM z punktu 1:

Ostrzeżenie : Na kod może być niechlujny/nieoptymalny, ponieważ nie jestem programistą Pythona. Ale robi swoje.

import numpy as np
import math

#### E-M Coin Toss Example as given in the EM tutorial paper by Do and Batzoglou* #### 

def get_mn_log_likelihood(obs,probs):
    """ Return the (log)likelihood of obs, given the probs"""
    # Multinomial Distribution Log PMF
    # ln (pdf)      =             multinomial coeff            *   product of probabilities
    # ln[f(x|n, p)] = [ln(n!) - (ln(x1!)+ln(x2!)+...+ln(xk!))] + [x1*ln(p1)+x2*ln(p2)+...+xk*ln(pk)]     

    multinomial_coeff_denom= 0
    prod_probs = 0
    for x in range(0,len(obs)): # loop through state counts in each observation
        multinomial_coeff_denom = multinomial_coeff_denom + math.log(math.factorial(obs[x]))
        prod_probs = prod_probs + obs[x]*math.log(probs[x])

    multinomial_coeff = math.log(math.factorial(sum(obs))) -  multinomial_coeff_denom
    likelihood = multinomial_coeff + prod_probs
    return likelihood

# 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((5,2), dtype=float) 
    expectation_B = np.zeros((5,2), dtype=float)
    for i in range(0,len(experiments)):
        e = experiments[i] # i'th experiment
        ll_A = get_mn_log_likelihood(e,np.array([pA_heads[j],1-pA_heads[j]])) # loglikelihood of e given coin A
        ll_B = get_mn_log_likelihood(e,np.array([pB_heads[j],1-pB_heads[j]])) # loglikelihood of e given coin B

        weightA = math.exp(ll_A) / ( math.exp(ll_A) + math.exp(ll_B) ) # corresponding weight of A proportional to likelihood of A 
        weightB = math.exp(ll_B) / ( math.exp(ll_A) + math.exp(ll_B) ) # corresponding weight of B proportional to likelihood of 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

Technicznie termin " EM " jest trochę niedopowiedziany, ale zakładam, że odnosisz się do techniki Gaussian Mixture Modeling cluster analysis, która jest instancją ogólnej zasady EM.

Właściwie, analiza klastra EM nie jest klasyfikatorem. Wiem, że niektórzy uważają klastry za "klasyfikację bez nadzoru", ale w rzeczywistości analiza klastrów jest czymś zupełnie innym.

Kluczowa różnica i wielkie nieporozumienie zawsze z analizy klastra jest to, że: w analaysis klastra, nie ma "poprawnego rozwiązania" . Jest to metoda wiedzy odkrycia , w rzeczywistości ma na celu znalezienie czegoś Nowego ! To sprawia, że ocena jest bardzo trudna. Często jest ona oceniana przy użyciu znanej klasyfikacji jako odniesienie, ale nie zawsze jest to właściwe: klasyfikacja, którą posiadasz, może lub nie odzwierciedlać tego, co znajduje się w danych.

Pozwól, że podam ci przykład: masz duży zestaw danych klientów, w tym dane dotyczące płci. Metoda, która dzieli ten zestaw danych na "męski" i "żeński" jest optymalna, jeśli porównasz go z istniejącymi klasami. W" przewidywania " sposób myślenia jest to dobre, jak dla nowych użytkowników można teraz przewidzieć ich płeć. W "odkrywaniu wiedzy" sposób myślenia jest naprawdę zły, ponieważ chciałeś odkryć jakąś nową strukturę w danych. Metoda, która np. podzieliłaby dane na osoby starsze i dzieci, jednak wynikałaby gorzej, jak to tylko możliwe zdobądź w odniesieniu do klasy męsko-żeńskiej. Byłby to jednak doskonały wynik grupowania (gdyby nie podano wieku).

A teraz wracamy do EM. Zasadniczo zakłada, że Twoje dane składają się z wielu wielowymiarowych rozkładów normalnych (zauważ, że jest to silne założenie Bardzo, w szczególności gdy ustalasz liczbę klastrów!). Następnie próbuje znaleźć lokalny model optymalny dla tego przez naprzemiennie ulepszając model i przypisanie obiektu do model .

Aby uzyskać najlepsze wyniki w kontekście klasyfikacji, Wybierz liczbę klastrów większych niż liczba klas, lub nawet zastosuj klastry tylko dopojedynczych klas (aby dowiedzieć się, czy istnieje jakaś struktura w klasie!).

Powiedz, że chcesz szkolić klasyfikatora, aby odróżnić "samochody", "rowery" i "ciężarówki". Przy założeniu, że dane składają się dokładnie z 3 rozkładów normalnych, nie ma większego zastosowania. Można jednak założyć, że jest więcej niż jeden rodzaj samochodów (oraz ciężarówek i rowerów). Zamiast więc szkolić klasyfikatora dla tych trzech klas, klastrujesz samochody, ciężarówki i rowery na 10 klastrów (a może 10 samochodów, 3 ciężarówki i 3 rowery, cokolwiek), a następnie szkolisz klasyfikatora, aby odróżnić te 30 klas, a następnie scalić wynik klasy z powrotem do oryginalnych klas. Można również odkryć, że istnieje jeden klaster, który jest szczególnie trudny do sklasyfikowania, na przykład Trikes. To trochę samochody, a trochę rowery. Lub Ciężarówki dostawcze, które są bardziej jak samochody ponadgabarytowe niż Ciężarówki.

Inne odpowiedzi będąc dobrym, postaram się przedstawić inną perspektywę i zająć się intuicyjną częścią pytania.

Algorytm EM (Expectation-Maximization) jest odmianą klasy algorytmów iteracyjnych wykorzystujących dualność

Fragment (podkreślenie):

W matematyce dualizm, ogólnie rzecz biorąc, tłumaczy pojęcia, twierdzeń lub struktur matematycznych do innych pojęć, twierdzeń lub struktur, w sposób indywidualny, często (ale nie zawsze) za pomocą operacji inwolucyjnej: jeśli dual Z A jest B, to dual Z B jest A. takie inwolucje czasami mają PUNKTY STAŁE , tak że podwójny of a Is A itself

Zazwyczaj a podwójny B obiektu A jest związany z a w pewien sposób, który zachowuje pewną symetrię lub zgodność . Na przykład AB = const

Przykłady algorytmów iteracyjnych, wykorzystujących dualność (w poprzednim znaczeniu) są:

1. algorytm Euklidesowy dla największego wspólnego dzielnika i jego wariantów
2. algorytm bazowy wektora grama-Schmidta i warianty
3. średnia arytmetyczna-średnia geometryczna nierówności i jej warianty
4. Algorytm maksymalizacji oczekiwań i jego warianty (Zobacz także tutaj, aby uzyskać widok geometryczny informacji )
5. (.. inne podobne algorytmy..)

W podobny sposób, EM algorytm może być również postrzegany jako dwa podwójne kroki maksymalizacji :

..[EM] jest postrzegana jako maksymalizacja funkcji wspólnych parametrów i rozkład na nieobserwowane zmienne.. E-step maksymalizuje funkcja ta w odniesieniu do rozkładu nad nieobserwowanym zmienne; M-krok w odniesieniu do parametrów..

W algorytmie iteracyjnym wykorzystującym dualność istnieje jednoznaczne (lub niejawne) założenie punktu równowagi (lub stałego) zbieżność (dla EM jest to udowodnione za pomocą nierówności Jensena)

Więc zarys takich algorytmów jest:

1. E-podobny krok: Znajdź najlepsze rozwiązanie x w odniesieniu do danej y trzymanej stałej.
2. M-jak krok (podwójny): Znajdź najlepsze rozwiązanie y w odniesieniu do x (jak obliczono w poprzednim kroku) jest utrzymywana stała.
3. kryterium zakończenia / stopnia zbieżności: powtórz kroki 1, 2 z zaktualizowane wartości x,y aż do osiągnięcia zbieżności (lub określonej liczby iteracji)

Zauważ, że gdy taki algorytm zbiega się do (globalnego) optimum, to znalazł konfigurację, która jest najlepsza w obu zmysłach (tj. zarówno w x domena/parametry i y domena / parametry). Jednak algorytm może po prostu znaleźć lokalne optimum, a nie globalne optimum.

Powiedziałbym jest to intuicyjny opis zarysu algorytmu

W przypadku argumentów statystycznych i zastosowań inne odpowiedzi dały dobre wyjaśnienia (sprawdź również odniesienia w tej odpowiedzi)

Przyjęta odpowiedź odwołuje się do Chuong em Paper , który porządnie tłumaczy em. Istnieje również wideo na youtube , które wyjaśnia artykuł bardziej szczegółowo.

Podsumowując, oto scenariusz:

1st:  {H,T,T,T,H,H,T,H,T,H} 5 Heads, 5 Tails; Did coin A or B generate me?
2nd:  {H,H,H,H,T,H,H,H,H,H} 9 Heads, 1 Tails
3rd:  {H,T,H,H,H,H,H,T,H,H} 8 Heads, 2 Tails
4th:  {H,T,H,T,T,T,H,H,T,T} 4 Heads, 6 Tails
5th:  {T,H,H,H,T,H,H,H,T,H} 7 Heads, 3 Tails

Two possible coins, A & B are used to generate these distributions.
A & B have an unknown parameter: their bias towards heads.

We don't know the biases, but we can simply start with a guess: A=60% heads, B=50% heads.

W przypadku pierwszego pytania próbnego, intuicyjnie sądzimy, że b je wygenerował, ponieważ proporcja głów bardzo dobrze pasuje do stronniczości B... ale ta wartość była tylko domysłem, więc nie możemy być pewni.

Mając to na uwadze, lubię myśleć o rozwiązanie EM takie jak to:

• każda próba rzutów "głosuje" na którą monetę lubi najbardziej
  • jest to oparte na tym, jak dobrze każda moneta pasuje do swojej dystrybucji
  • lub, z punktu widzenia monety, istnieje wysokie oczekiwanie na zobaczenie tej próby w stosunku do drugiej monety (na podstawielog likelihoods ).
• w zależności od tego, jak bardzo każda próba lubi każdą monetę, może zaktualizować domyślanie parametru danej monety (bias).
  • im bardziej proces lubi monetę,tym bardziej aktualizuje stronniczość monety, aby odzwierciedlić jej własne!
  • zasadniczo błędy monety są aktualizowane przez połączenie tych ważonych aktualizacji we wszystkich próbach, proces zwany (maximazation {20]}), który odnosi się do próby uzyskania najlepszych domysłów dla błędów każdej monety, biorąc pod uwagę zestaw prób.

EM jest używany do maksymalizacji prawdopodobieństwa modelu Q ze zmiennymi utajonymi Z.

To iteracyjna optymalizacja.

theta <- initial guess for hidden parameters
while not converged:
    #e-step
    Q(theta'|theta) = E[log L(theta|Z)]
    #m-step
    theta <- argmax_theta' Q(theta'|theta)

E-krok: biorąc pod uwagę bieżące oszacowanie Z Oblicz oczekiwaną funkcję loglikelihood

M-krok: znajdź theta, która maksymalizuje to Q

Przykład GMM:

E-step: szacowanie przydziałów etykiet dla każdego punktu danych, biorąc pod uwagę aktualną estymację parametru gmm

M-step: assigments

K-means jest również algorytmem EM i istnieje wiele wyjaśniających animacji na k-means.

Używając tego samego artykułu Do I Batzoglou cytowanego w odpowiedzi Zhubarba, zaimplementowałem EM dla tego problemu wJava . Komentarze do jego odpowiedzi pokazują, że algorytm utkwi w lokalnym optimum, co występuje również w mojej implementacji, jeśli parametry thetaA i thetaB są takie same.

Poniżej znajduje się standardowe wyjście mojego kodu, pokazujące zbieżność parametrów.

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

Poniżej moja Java implementacja EM do rozwiązania problemu w (Do I Batzoglou, 2008). Podstawową częścią implementacji jest pętla do uruchamiania EM aż do zbieżności parametrów.

private Parameters _parameters;

public Parameters run()
{
    while (true)
    {
        expectation();

        Parameters estimatedParameters = maximization();

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

        _parameters = estimatedParameters;
    }

    return _parameters;
}

Poniżej znajduje się cały kod.

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,
respectively.
*****************************************************************************/
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
experimental 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);
    }
}