La maximisation des attentes (EM) est une sorte de méthode probabiliste pour classer les données. Veuillez me corriger si je me trompe si ce n'est pas un classificateur.

Qu'est-ce qu'une explication intuitive de cette technique EM? Qu'y a-t-il expectationici et qu'est-ce que l'être maximized?

Remarque: le code derrière cette réponse se trouve ici .

Supposons que nous ayons des données échantillonnées dans deux groupes différents, rouge et bleu:

Ici, nous pouvons voir quel point de données appartient au groupe rouge ou bleu. Cela facilite la recherche des paramètres qui caractérisent chaque groupe. Par exemple, la moyenne du groupe rouge est d'environ 3, la moyenne du groupe bleu est d'environ 7 (et nous pourrions trouver les moyennes exactes si nous le voulions).

C'est ce qu'on appelle généralement l' estimation du maximum de vraisemblance . Compte tenu de certaines données, nous calculons la valeur d'un paramètre (ou de paramètres) qui explique le mieux ces données.

Imaginez maintenant que nous ne pouvons pas voir quelle valeur a été échantillonnée dans quel groupe. Tout nous semble violet:

Ici, nous savons qu'il existe deux groupes de valeurs, mais nous ne savons pas à quel groupe appartient une valeur particulière.

Pouvons-nous encore estimer les moyennes du groupe rouge et du groupe bleu qui correspondent le mieux à ces données?

Oui, souvent nous le pouvons! La maximisation des attentes nous donne un moyen de le faire. L'idée très générale derrière l'algorithme est la suivante:

1. Commencez par une première estimation de ce que pourrait être chaque paramètre.
2. Calculez la probabilité que chaque paramètre produise le point de données.
3. Calculez les poids pour chaque point de données en indiquant s'il est plus rouge ou plus bleu en fonction de la probabilité qu'il soit produit par un paramètre. Combinez les pondérations avec les données ( attente ).
4. Calculez une meilleure estimation des paramètres en utilisant les données ajustées en poids ( maximisation ).
5. Répétez les étapes 2 à 4 jusqu'à ce que l'estimation du paramètre converge (le processus arrête de produire une estimation différente).

Ces étapes nécessitent des explications supplémentaires, je vais donc parcourir le problème décrit ci-dessus.

Exemple: estimation de la moyenne et de l'écart type

J'utiliserai Python dans cet exemple, mais le code devrait être assez facile à comprendre si vous n'êtes pas familier avec ce langage.

Supposons que nous ayons deux groupes, rouge et bleu, avec les valeurs distribuées comme dans l'image ci-dessus. Plus précisément, chaque groupe contient une valeur tirée d'une distribution normale avec les paramètres suivants:

import numpy as np
from scipy import stats

np.random.seed(110) # for reproducible 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))) # for later use...

Voici à nouveau une image de ces groupes rouges et bleus (pour vous éviter d'avoir à faire défiler vers le haut):

Lorsque nous pouvons voir la couleur de chaque point (c'est-à-dire à quel groupe il appartient), il est très facile d'estimer la moyenne et l'écart type de chaque groupe. Nous transmettons simplement les valeurs rouge et bleu aux fonctions intégrées de NumPy. Par exemple:

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

Mais que faire si nous ne pouvons pas voir les couleurs des points? Autrement dit, au lieu de rouge ou de bleu, chaque point a été coloré en violet.

Pour essayer de récupérer les paramètres d'écart moyen et standard pour les groupes rouge et bleu, nous pouvons utiliser la maximisation des attentes.

Notre première étape ( étape 1 ci-dessus) consiste à deviner les valeurs des paramètres pour la moyenne et l'écart type de chaque groupe. Nous n'avons pas à deviner intelligemment; nous pouvons choisir tous les nombres que nous aimons:

# 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

Ces estimations de paramètres produisent des courbes en cloche qui ressemblent à ceci:

Ce sont de mauvaises estimations. Les deux moyens (les lignes pointillées verticales) semblent loin de tout type de «milieu» pour des groupes de points sensibles, par exemple. Nous voulons améliorer ces estimations.

L'étape suivante ( étape 2 ) consiste à calculer la probabilité que chaque point de données apparaisse sous le paramètre actuel suppose:

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)

Ici, nous avons simplement mis chaque point de données dans la fonction de densité de probabilité pour une distribution normale en utilisant nos estimations actuelles à la moyenne et à l'écart type pour le rouge et le bleu. Cela nous indique, par exemple, qu'avec nos estimations actuelles, le point de données à 1,761 est beaucoup plus susceptible d'être rouge (0,189) que bleu (0,00003).

Pour chaque point de données, nous pouvons transformer ces deux valeurs de vraisemblance en poids ( é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 pour la moyenne et l'écart type des groupes rouge et bleu ( étape 4 ).

Nous calculons deux fois la moyenne et l'écart type en utilisant tous les points de données, mais avec des pondérations différentes: une fois pour les poids rouges et une fois pour les poids bleus.

L'intuition clé est que plus le poids d'une couleur sur un point de données est élevé, plus le point de données influence les estimations suivantes pour les paramètres de cette couleur. Cela a pour effet de "tirer" les paramètres dans la bonne direction.

def estimate_mean(data, weight):
    """
    For each data point, multiply the point by the probability it
    was drawn from the colour's distribution (its "weight").

    Divide by the total weight: essentially, we're finding where 
    the weight is centred among our data points.
    """
    return np.sum(data * weight) / np.sum(weight)

def estimate_std(data, weight, mean):
    """
    For each data point, multiply the point's squared difference
    from a mean value by the probability it was drawn from
    that distribution (its "weight").

    Divide by the total weight: essentially, we're finding where 
    the weight is centred among the values for the difference of
    each data point from the mean.

    This is the estimate of the variance, take the positive square
    root to find the standard deviation.
    """
    variance = np.sum(weight * (data - mean)**2) / np.sum(weight)
    return np.sqrt(variance)

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

Nous avons de nouvelles estimations pour les paramètres. Pour les améliorer à nouveau, nous pouvons revenir à l'étape 2 et répéter le processus. Nous faisons cela jusqu'à ce que les estimations convergent, ou après qu'un certain nombre d'itérations aient été effectuées ( étape 5 ).

Pour nos données, les cinq premières itérations de ce processus ressemblent à ceci (les itérations récentes ont une apparence plus forte):

On voit que les moyennes convergent déjà sur certaines valeurs, et les formes des courbes (régies par l'écart type) deviennent également plus stables.

Si nous continuons pendant 20 itérations, nous nous retrouvons avec ce qui suit:

Le processus EM a convergé vers les valeurs suivantes, qui s'avèrent très proches des valeurs réelles (où nous pouvons voir les couleurs - pas de variables cachées):

          | EM guess | Actual |  Delta
----------+----------+--------+-------
Red mean  |    2.910 |  2.802 |  0.108
Red std   |    0.854 |  0.871 | -0.017
Blue mean |    6.838 |  6.932 | -0.094
Blue std  |    2.227 |  2.195 |  0.032

Dans le code ci-dessus, vous avez peut-être remarqué que la nouvelle estimation de l'écart type a été calculée en utilisant l'estimation de l'itération précédente pour la moyenne. En fin de compte, peu importe si nous calculons d'abord une nouvelle valeur pour la moyenne, car nous cherchons simplement la variance (pondérée) des valeurs autour d'un point central. Nous verrons toujours les estimations des paramètres converger.

EM est un algorithme pour maximiser une fonction de vraisemblance lorsque certaines des variables de votre modèle ne sont pas observées (c'est-à-dire lorsque vous avez des variables latentes).

Vous pourriez vous demander, si nous essayons simplement de maximiser une fonction, pourquoi n'utilisons-nous pas simplement la machinerie existante pour maximiser une fonction. Eh bien, si vous essayez de maximiser cela en prenant des dérivées et en les mettant à zéro, vous constatez que dans de nombreux cas, les conditions du premier ordre n'ont pas de solution. Il y a un problème de poule et d'œuf en ce sens que pour résoudre les paramètres de votre modèle, vous devez connaître la distribution de vos données non observées; mais la distribution de vos données non observées est fonction des paramètres de votre modèle.

EM essaie de contourner ce problème en devinant itérativement une distribution pour les données non observées, puis en estimant les paramètres du modèle en maximisant quelque chose qui est une limite inférieure de la fonction de vraisemblance réelle, et en répétant jusqu'à la convergence:

L'algorithme EM

Commencez par deviner les valeurs de vos paramètres de modèle

Étape E: Pour chaque point de données qui a des valeurs manquantes, utilisez votre équation de modèle pour résoudre la distribution des données manquantes compte tenu de votre estimation actuelle des paramètres du modèle et compte tenu des données observées (notez que vous résolvez une distribution pour chaque valeur, pas pour la valeur attendue). Maintenant que nous avons une distribution pour chaque valeur manquante, nous pouvons calculer l' espérance de la fonction de vraisemblance par rapport aux variables non observées. Si notre estimation du paramètre du modèle était correcte, cette probabilité attendue sera la vraisemblance de nos données observées; si les paramètres n'étaient pas corrects, ce ne sera qu'une borne inférieure.

Étape M: Maintenant que nous avons une fonction de vraisemblance attendue sans variables non observées, maximisez la fonction comme vous le feriez dans le cas pleinement observé, pour obtenir une nouvelle estimation des paramètres de votre modèle.

Répétez jusqu'à convergence.

Voici une recette simple pour comprendre l'algorithme de maximisation des attentes:

1- Lisez ce tutoriel EM par Do et Batzoglou.

2- Vous pouvez avoir des points d'interrogation dans votre tête, jetez un œil aux explications sur cette page d' échange de pile de maths .

3- Regardez ce code que j'ai écrit en Python qui explique l'exemple du tutoriel EM de l'élément 1:

Attention: le code peut être désordonné / sous-optimal, car je ne suis pas un développeur Python. Mais ça fait le travail.

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

Techniquement, le terme «EM» est un peu sous-spécifié, mais je suppose que vous faites référence à la technique d'analyse en grappes de modélisation de mélange gaussien, qui est un exemple du principe général EM.

En fait, l' analyse de cluster EM n'est pas un classificateur . Je sais que certaines personnes considèrent le regroupement comme une «classification non supervisée», mais en fait, l'analyse de cluster est quelque chose de tout à fait différent.

La principale différence, et le grand malentendu que les gens de classification ont toujours avec l'analyse par grappes, est que: dans l' analyse par grappes, il n'y a pas de «solution correcte» . C'est une méthode de découverte de connaissances , elle vise en fait à trouver quelque chose de nouveau ! Cela rend l'évaluation très délicate. Il est souvent évalué en utilisant une classification connue comme référence, mais ce n'est pas toujours approprié: la classification que vous avez peut refléter ou non ce que contiennent les données.

Permettez-moi de vous donner un exemple: vous avez un grand ensemble de données de clients, y compris des données de genre. Une méthode qui divise cet ensemble de données en «homme» et «femme» est optimale lorsque vous la comparez avec les classes existantes. Dans une façon de penser «prédictive», c'est bien, car pour les nouveaux utilisateurs, vous pouvez désormais prédire leur sexe. Dans une manière de penser «découverte de connaissances», c'est en fait mauvais, car vous vouliez découvrir une nouvelle structure dans les données. Une méthode qui, par exemple, diviserait les données en personnes âgées et en enfants, obtiendrait cependant un score aussi pire que possible par rapport à la classe homme / femme. Cependant, ce serait un excellent résultat de regroupement (si l'âge n'était pas indiqué).

Revenons maintenant à EM. Essentiellement, cela suppose que vos données sont composées de plusieurs distributions normales multivariées (notez qu'il s'agit d'une hypothèse très forte, en particulier lorsque vous fixez le nombre de clusters!). Il essaie ensuite de trouver un modèle local optimal pour cela en améliorant alternativement le modèle et l'affectation des objets au modèle .

Pour de meilleurs résultats dans un contexte de classification, choisissez le nombre de clusters supérieur au nombre de classes, ou même appliquez le clustering à des classes uniques uniquement (pour savoir s'il existe une structure au sein de la classe!).

Supposons que vous souhaitiez former un classificateur à distinguer les «voitures», les «vélos» et les «camions». Il est peu utile de supposer que les données consistent exactement en 3 distributions normales. Cependant, vous pouvez supposer qu'il existe plus d'un type de voitures (et camions et vélos). Donc, au lieu de former un classificateur pour ces trois classes, vous regroupez les voitures, les camions et les vélos en 10 groupes chacun (ou peut-être 10 voitures, 3 camions et 3 vélos, peu importe), puis formez un classificateur pour distinguer ces 30 classes, puis fusionne le résultat de la classe avec les classes d'origine. Vous pouvez également découvrir qu'il existe un cluster qui est particulièrement difficile à classer, par exemple Trikes. Ce sont un peu des voitures et un peu des vélos. Ou des camions de livraison, qui ressemblent plus à des voitures surdimensionnées qu'à des camions.

Les autres réponses étant bonnes, je vais essayer de fournir une autre perspective et d'aborder la partie intuitive de la question.

L'algorithme EM (Expectation-Maximization) est une variante d'une classe d'algorithmes itératifs utilisant la dualité

Extrait (c'est moi qui souligne):

En mathématiques, une dualité, d'une manière générale, traduit des concepts, théorèmes ou structures mathématiques en d'autres concepts, théorèmes ou structures, de manière biunivoque, souvent (mais pas toujours) au moyen d'une opération d'involution: si le duel de A est B, alors le dual de B est A. De telles involutions ont parfois des points fixes , de sorte que le dual de A est A lui-même

Habituellement, un double B d'un objet A est lié à A d'une manière qui préserve une certaine symétrie ou compatibilité . Par exemple AB = const

Des exemples d'algorithmes itératifs, employant la dualité (au sens précédent) sont:

1. Algorithme euclidien pour le plus grand diviseur commun et ses variantes
2. Algorithme de base vectorielle de Gram – Schmidt et variantes
3. Moyenne arithmétique - Inégalité de la moyenne géométrique et ses variantes
4. Algorithme de maximisation des attentes et ses variantes (voir aussi ici pour une vue géométrique d'information )
5. (.. autres algorithmes similaires ..)

De la même manière, l'algorithme EM peut également être vu comme deux étapes de maximisation doubles :

.. [EM] est considéré comme maximisant une fonction conjointe des paramètres et de la distribution sur les variables non observées. L'étape E maximise cette fonction par rapport à la distribution sur les variables non observées; le pas M par rapport aux paramètres.

Dans un algorithme itératif utilisant la dualité, il y a l'hypothèse explicite (ou implicite) d'un point de convergence d'équilibre (ou fixe) (pour EM, cela est prouvé en utilisant l'inégalité de Jensen)

Ainsi, le contour de ces algorithmes est:

1. Étape de type E: Trouver la meilleure solution x par rapport à y donné étant maintenu constant.
2. Étape de type M (double): Trouver la meilleure solution y par rapport à x (tel que calculé à l'étape précédente) maintenu constant.
3. Critère d'étape de terminaison / convergence: Répétez les étapes 1, 2 avec les valeurs mises à jour de x , y jusqu'à ce que la convergence (ou le nombre spécifié d'itérations soit atteint)

Notez que lorsqu'un tel algorithme converge à un (global) optimal, il a trouvé une configuration qui est mieux dans les deux sens (dans les deux x domaine / paramètres et les y domaine / paramètres). Cependant, l'algorithme peut simplement trouver un optimum local et non l' optimum global .

je dirais que c'est la description intuitive du contour de l'algorithme

Pour les arguments statistiques et les applications, d'autres réponses ont donné de bonnes explications (vérifiez également les références dans cette réponse)

La réponse acceptée fait référence au Chuong EM Paper , qui fait un travail décent pour expliquer la SE. Il existe également une vidéo sur youtube qui explique l'article plus en détail.

Pour récapituler, voici le scénario:

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.

Dans le cas de la question du premier essai, intuitivement, nous pensons que B l'a générée car la proportion de têtes correspond très bien au biais de B ... mais cette valeur n'était qu'une supposition, nous ne pouvons donc pas en être sûrs.

Dans cet esprit, j'aime penser à la solution EM comme ceci:

• Chaque essai de flips permet de `` voter '' sur la pièce qui lui plaît le plus
  • Ceci est basé sur l'adéquation de chaque pièce à sa distribution
  • OU, du point de vue de la pièce, il y a de fortes attentes de voir ce procès par rapport à l'autre pièce (sur la base des probabilités logarithmiques ).
• En fonction de combien chaque essai aime chaque pièce, il peut mettre à jour l'estimation du paramètre de cette pièce (biais).
  • Plus un essai aime une pièce, plus il met à jour le biais de la pièce pour refléter le sien!
  • Essentiellement, les biais de la pièce sont mis à jour en combinant ces mises à jour pondérées dans tous les essais, un processus appelé ( maximazation ), qui consiste à essayer d'obtenir les meilleures estimations pour le biais de chaque pièce compte tenu d'un ensemble d'essais.

Cela peut être une simplification excessive (ou même fondamentalement faux à certains niveaux), mais j'espère que cela aidera à un niveau intuitif!

EM est utilisé pour maximiser la probabilité d'un modèle Q avec des variables latentes Z.

C'est une optimisation itérative.

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-step: étant donné l'estimation actuelle de Z, calculer la fonction log-vraisemblance attendue

m-step: trouver thêta qui maximise ce Q

Exemple GMM:

e-step: estimer les attributions d'étiquettes pour chaque point de données en fonction de l'estimation actuelle du paramètre gmm

m-step: maximiser un nouveau thêta compte tenu des nouvelles attributions d'étiquettes

K-means est également un algorithme EM et il y a beaucoup d'animations explicatives sur K-means.

En utilisant le même article de Do et Batzoglou cité dans la réponse de Zhubarb, j'ai implémenté EM pour ce problème en Java . Les commentaires à sa réponse montrent que l'algorithme se bloque à un optimum local, ce qui se produit également avec mon implémentation si les paramètres thetaA et thetaB sont les mêmes.

Voici la sortie standard de mon code, montrant la convergence des paramètres.

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

Voici mon implémentation Java de EM pour résoudre le problème dans (Do et Batzoglou, 2008). La partie centrale de l'implémentation est la boucle pour exécuter EM jusqu'à ce que les paramètres convergent.

private Parameters _parameters;

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

        Parameters estimatedParameters = maximization();

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

        _parameters = estimatedParameters;
    }

    return _parameters;
}

Voici le code complet.

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