Exemple numérique pour comprendre l'attente-maximisation

J'essaie de bien comprendre l'algorithme EM, de pouvoir l'implémenter et de l'utiliser. J'ai passé une journée entière à lire la théorie et un article dans lesquels EM est utilisé pour suivre un avion en utilisant les informations de position provenant d'un radar. Honnêtement, je ne pense pas comprendre parfaitement l’idée sous-jacente. Quelqu'un peut-il m'indiquer un exemple numérique montrant quelques itérations (3-4) de l'EM pour un problème plus simple (comme l'estimation des paramètres d'une distribution gaussienne ou d'une séquence d'une série sinusoïdale ou l'ajustement d'une ligne).

Même si quelqu'un peut m'indiquer un morceau de code (avec des données synthétiques), je peux essayer de le parcourir.

Voici une recette pour apprendre la SE avec un exemple pratique et (à mon avis) très intuitif:

1. Lisez ce court didacticiel sur l’EM de Do et Batzoglou. C'est le schéma où l'exemple de tirage au sort est expliqué:

   

2. Vous avez peut-être des points d'interrogation dans la tête, notamment en ce qui concerne l'origine des probabilités de l'étape des attentes. Veuillez consulter les explications de cette page d’ échange de piles mathématiques .

3. Regardez / exécutez ce code que j'ai écrit en Python et qui simule la solution au problème du tirage au sort dans l'article du didacticiel sur la ME du point 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()

Il semble que votre question comporte deux parties: l’idée sous-jacente et un exemple concret. Je vais commencer par l'idée sous-jacente, puis un lien vers un exemple en bas.

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

Le cas le plus courant auquel les gens s’occupent est probablement les distributions mixtes. Pour notre exemple, examinons un modèle de mélange gaussien simple:

Vous avez deux distributions gaussiennes univariées différentes avec des moyennes et des variances d'unités différentes.

Vous avez toute une série de points de données, mais vous ne savez pas quels points proviennent de quelle distribution et vous ne savez pas non plus quelle est la moyenne des deux distributions.

Et maintenant tu es coincé:

• Si vous connaissiez le vrai moyen, vous pourriez déterminer quels points de données provenaient de quel gaussien. Par exemple, si un point de données a une valeur très élevée, il provient probablement de la distribution avec la moyenne la plus élevée. Mais vous ne savez pas quels sont les moyens, alors cela ne fonctionnera pas.

• Si vous saviez de quelle distribution provenait chaque point, vous pourriez alors estimer la moyenne des deux distributions à l'aide de la moyenne d'échantillonnage des points pertinents. Mais vous ne savez pas réellement quels points attribuer à quelle distribution, cela ne fonctionnera donc pas non plus.

Donc, aucune de ces approches ne semble fonctionner: vous auriez besoin de connaître la réponse avant de pouvoir trouver la réponse, et vous êtes bloqué.

Ce que vous permet de faire est d’alterner entre ces deux étapes faciles au lieu de s’occuper de tout le processus à la fois.

Vous devrez commencer par deviner les deux moyens (bien que votre estimation ne doive pas nécessairement être très précise, vous devez commencer quelque part).

Si votre hypothèse concernant les moyennes était exacte, vous auriez alors suffisamment d'informations pour exécuter l'étape de mon premier point ci-dessus, et vous pourriez (probablement) attribuer chaque point de données à l'un des deux Gaussiens. Même si nous savons que notre hypothèse est fausse, essayons quand même. Et puis, étant donné les distributions attribuées à chaque point, vous pouvez obtenir de nouvelles estimations pour les moyennes en utilisant le deuxième point. Il s'avère que chaque fois que vous parcourez ces deux étapes, vous améliorez la limite inférieure de la probabilité du modèle.

C'est déjà assez cool: même si les deux suggestions dans les points ci-dessus ne semblaient pas fonctionner individuellement, vous pouvez toujours les utiliser ensemble pour améliorer le modèle. La vraie magie de l'EM est que, après suffisamment d'itérations, la limite inférieure sera si élevée qu'il n'y aura plus d'espace entre elle et le maximum local. En conséquence, et vous avez optimisé localement la probabilité.

Ainsi, vous n'avez pas seulement amélioré le modèle, vous avez trouvé le meilleur modèle possible avec les mises à jour incrémentielles.

Cette page de Wikipedia montre un exemple un peu plus compliqué (Gaussiennes à deux dimensions et covariance inconnue), mais l’idée de base est la même. Il inclut également du Rcode bien commenté pour implémenter l'exemple.

Dans le code, l'étape "Expectation" (étape E) correspond à mon premier point: déterminer quel gaussien assume la responsabilité de chaque point de données, compte tenu des paramètres actuels de chaque gaussien. L'étape "Maximisation" (étape M) met à jour les moyennes et les covariances, compte tenu de ces affectations, comme dans mon deuxième point.

Comme vous pouvez le constater dans l'animation, ces mises à jour permettent rapidement à l'algorithme de passer d'un ensemble d'estimations terribles à un ensemble d'excellentes: il semble vraiment y avoir deux nuages ​​de points centrés sur les deux distributions gaussiennes trouvées par EM.

Voici un exemple d'optimisation des attentes (EM) utilisé pour estimer la moyenne et l'écart type. Le code est en Python, mais il devrait être facile à suivre même si vous n'êtes pas familier avec le langage.

La motivation pour l'EM

Les points rouges et bleus ci-dessous sont tirés de deux distributions normales différentes, chacune avec une moyenne et un écart type particuliers:

Pour calculer des approximations raisonnables des paramètres "vrais" de moyenne et d'écart type pour la distribution rouge, nous pourrions très facilement examiner les points rouges et enregistrer la position de chacun d'eux, puis utiliser les formules familières (et de la même manière pour le groupe bleu). .

Considérons maintenant le cas où nous savons qu’il existe deux groupes de points, mais nous ne pouvons pas voir quel point appartient à quel groupe. En d'autres termes, les couleurs sont cachées:

Il n’est pas du tout évident de diviser les points en deux groupes. Nous ne sommes plus en mesure de regarder les positions et de calculer des estimations pour les paramètres de la distribution rouge ou de la distribution bleue.

C’est là que l’EM peut être utilisé pour résoudre le problème.

Utiliser EM pour estimer les paramètres

Voici le code utilisé pour générer les points indiqués ci-dessus. Vous pouvez voir les moyennes et les écarts-types réels des distributions normales à partir desquelles les points ont été tirés. Les variables redet bluetiennent les positions de chaque point dans les groupes rouge et bleu respectivement:

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

Si nous pouvions voir la couleur de chaque point, nous essaierions de récupérer les moyennes et les écarts-types à l'aide des fonctions de la bibliothèque:

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

Mais puisque les couleurs nous sont cachées, nous allons commencer le processus EM ...

Premièrement, nous devinons simplement les valeurs des paramètres de chaque groupe ( étape 1 ). Ces suppositions ne doivent pas forcément être bonnes:

# 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

Assez mauvaises conjectures - les moyens semblent être loin du "milieu" d'un groupe de points.

Pour continuer avec EM et améliorer ces suppositions, nous calculons la probabilité que chaque point de données (quelle que soit sa couleur secrète) apparaisse sous ces suppositions pour la moyenne et l'écart type ( étape 2 ).

La variable both_colourscontient chaque point de données. La fonction stats.normcalcule la probabilité du point sous une distribution normale avec les paramètres donnés:

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)

Cela nous indique, par exemple, qu’avec nos suppositions actuelles, le point de données situé à 1,761 est beaucoup plus susceptible d’être rouge (0,189) que bleu (0,00003).

Nous pouvons transformer ces deux valeurs de vraisemblance en pondérations ( étape 3 ) afin qu'elles totalisent 1 comme suit:

likelihood_total = likelihood_of_red + likelihood_of_blue

red_weight = likelihood_of_red / likelihood_total
blue_weight = likelihood_of_blue / likelihood_total

Avec nos estimations actuelles et nos poids nouvellement calculés, nous pouvons maintenant calculer de nouvelles estimations, probablement meilleures, pour les paramètres ( étape 4 ). Nous avons besoin d’une fonction pour la moyenne et d’une fonction pour l’écart type:

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)

Celles-ci ressemblent beaucoup aux fonctions habituelles à la moyenne et à l'écart type des données. La différence est l'utilisation d'unweight paramètre qui attribue une pondération à chaque point de données.

Cette pondération est la clé de la SE. Plus le poids d'une couleur sur un point de données est important, plus le point de données influence les estimations suivantes pour les paramètres de cette couleur. En fin de compte, cela a pour effet de tirer chaque paramètre dans la bonne direction.

Les nouvelles suppositions sont calculées avec ces fonctions:

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

Le processus EM est ensuite répété avec ces nouvelles suppositions à partir de l'étape 2. Nous pouvons répéter les étapes pour un nombre donné d'itérations (par exemple 20) ou jusqu'à ce que les paramètres convergent.

Après cinq itérations, nous voyons que nos mauvaises hypothèses initiales commencent à s’améliorer:

Après 20 itérations, le processus EM a plus ou moins convergé:

À des fins de comparaison, voici les résultats du processus EM comparés aux valeurs calculées pour lesquelles les informations de couleur ne sont pas masquées:

          | 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

Remarque: cette réponse a été adaptée de ma réponse sur Stack Overflow ici .

Suite à la réponse de Zhubarb, j'ai implémenté l'exemple EM "Lancer de la monnaie" de Do et Batzoglou dans GNU R. Notez que j'utilise la mlefonction du stats4package - cela m'a aidé à comprendre plus clairement le lien entre EM et 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))

Tout ce qui précède semble être une excellente ressource, mais je dois faire un lien avec ce bel exemple. Il présente une explication très simple pour trouver les paramètres de deux lignes d’un ensemble de points. Le tutoriel est de Yair Weiss au MIT.

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

La réponse donnée par Zhubarb est excellente, mais malheureusement, elle l’est en python. Ci-dessous, une implémentation Java de l'algorithme EM exécuté sur le même problème (posée dans l'article de Do et Batzoglou, 2008). J'ai ajouté quelques printf à la sortie standard pour voir comment les paramètres convergent.

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

Le code Java suit ci-dessous:

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

Eh bien, je vous suggérerais de lire un livre sur R de Maria L. Rizzo. L'un des chapitres contient l'utilisation de l'algorithme EM avec un exemple numérique. Je me souviens d’avoir parcouru le code pour une meilleure compréhension.

Aussi, essayez de le voir du point de vue du clustering au début. Élaborez à la main un problème de regroupement où 10 observations sont prises à partir de deux densités normales différentes. Cela devrait aider. Prenez l’aide de 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