Contoh numerik untuk memahami Ekspektasi-Maksimalisasi

Saya mencoba untuk mendapatkan pemahaman yang baik tentang algoritma EM, untuk dapat mengimplementasikan dan menggunakannya. Saya menghabiskan satu hari penuh membaca teori dan kertas di mana EM digunakan untuk melacak pesawat menggunakan informasi posisi yang berasal dari radar. Sejujurnya, saya pikir saya tidak sepenuhnya memahami ide yang mendasarinya. Adakah yang bisa mengarahkan saya ke contoh numerik yang menunjukkan beberapa iterasi (3-4) EM untuk masalah yang lebih sederhana (seperti memperkirakan parameter distribusi Gaussian atau urutan seri sinusoidal atau pemasangan garis).

Bahkan jika seseorang dapat mengarahkan saya ke sepotong kode (dengan data sintetis), saya dapat mencoba untuk menelusuri kode.

Ini adalah resep untuk belajar EM dengan contoh 'Coin-Toss' yang praktis dan (menurut saya) sangat intuitif:

1. Baca makalah tutorial EM singkat ini oleh Do dan Batzoglou. Ini adalah skema tempat contoh lemparan koin dijelaskan:

   

2. Anda mungkin memiliki tanda tanya di kepala Anda, terutama mengenai dari mana probabilitas dalam langkah Ekspektasi berasal. Silakan lihat penjelasan di halaman pertukaran tumpukan matematika ini .

3. Lihat / jalankan kode ini yang saya tulis dengan Python yang mensimulasikan solusi untuk masalah coin-toss di kertas tutorial EM item 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()
   

Sepertinya pertanyaan Anda memiliki dua bagian: ide yang mendasarinya dan contoh konkret. Saya akan mulai dengan ide yang mendasarinya, kemudian tautan ke contoh di bagian bawah.

EM berguna dalam Catch-22 situasi di mana sepertinya Anda perlu tahu sebelum Anda dapat menghitung dan Anda perlu tahu sebelum Anda dapat menghitung .B B $A$$B$$B$$A$

Kasus paling umum yang dihadapi adalah distribusi campuran. Sebagai contoh, mari kita lihat model campuran Gaussian sederhana:

Anda memiliki dua distribusi Gaussian univariat yang berbeda dengan cara dan varian unit yang berbeda.

Anda memiliki banyak titik data, tetapi Anda tidak yakin titik mana yang berasal dari distribusi mana, dan Anda juga tidak yakin tentang cara kedua distribusi.

Dan sekarang Anda terjebak:

• Jika Anda tahu cara yang sebenarnya, Anda bisa mengetahui titik data mana yang berasal dari mana Gaussian. Misalnya, jika suatu titik data memiliki nilai yang sangat tinggi, itu mungkin berasal dari distribusi dengan rata-rata yang lebih tinggi. Tetapi Anda tidak tahu apa artinya itu, jadi ini tidak akan berhasil.

• Jika Anda tahu dari mana asal masing-masing titik distribusi, maka Anda dapat memperkirakan sarana dua distribusi dengan menggunakan alat sampel dari titik-titik yang relevan. Tetapi Anda tidak benar-benar tahu poin mana yang harus ditetapkan untuk distribusi mana, jadi ini juga tidak akan berhasil.

Jadi tidak ada pendekatan yang berhasil: Anda harus tahu jawabannya sebelum Anda bisa menemukan jawabannya, dan Anda buntu.

Apa yang memungkinkan Anda lakukan EM adalah bergantian antara dua langkah yang bisa dilakukan ini alih-alih menangani seluruh proses sekaligus.

Anda harus mulai dengan menebak tentang dua cara (meskipun tebakan Anda tidak harus sangat akurat, Anda harus memulai suatu tempat).

Jika tebakan Anda tentang cara-cara itu akurat, maka Anda akan memiliki cukup informasi untuk melakukan langkah di poin pertama saya di atas, dan Anda dapat (secara probabilistik) menetapkan setiap titik data ke salah satu dari dua Gaussians. Meskipun kita tahu dugaan kita salah, mari kita coba saja. Dan kemudian, mengingat distribusi masing-masing titik yang ditentukan, Anda bisa mendapatkan perkiraan baru untuk sarana menggunakan poin kedua. Ternyata, setiap kali Anda melakukan loop melalui dua langkah ini, Anda meningkatkan batas bawah pada kemungkinan model.

Itu sudah sangat keren: meskipun kedua saran di poin-poin di atas sepertinya tidak akan bekerja secara individual, Anda masih dapat menggunakannya bersama-sama untuk meningkatkan model. The nyata keajaiban EM adalah bahwa, setelah cukup iterasi, batas bawah akan sangat tinggi sehingga ada tidak akan ada ruang antara itu dan maksimum lokal. Akibatnya, dan Anda telah mengoptimalkan kemungkinan secara lokal.

Jadi Anda belum memperbaiki model, Anda telah menemukan model terbaik yang dapat ditemukan dengan pembaruan tambahan.

Halaman ini dari Wikipedia menunjukkan contoh yang sedikit lebih rumit (Gaussians dua dimensi dan kovarian yang tidak diketahui), tetapi ide dasarnya sama. Ini juga termasuk Rkode yang dikomentari dengan baik untuk menerapkan contoh.

Dalam kode, langkah "Ekspektasi" (E-langkah) sesuai dengan poin-poin pertama saya: mencari tahu Gaussian mana yang bertanggung jawab untuk setiap titik data, mengingat parameter saat ini untuk setiap Gaussian. Langkah "Maksimalisasi" (langkah-M) memperbarui sarana dan kovarian, yang diberikan tugas ini, seperti pada poin kedua saya.

Seperti yang dapat Anda lihat dalam animasi, pembaruan ini dengan cepat memungkinkan algoritme beralih dari sekumpulan taksiran mengerikan ke sekumpulan yang sangat bagus: tampaknya memang ada dua titik awan yang berpusat pada dua distribusi Gaussian yang ditemukan EM.

Berikut adalah contoh dari Maksimalisasi Ekspektasi (EM) yang digunakan untuk memperkirakan mean dan standar deviasi. Kode ini dalam Python, tetapi harus mudah diikuti bahkan jika Anda tidak terbiasa dengan bahasa tersebut.

Motivasi untuk EM

Titik merah dan biru yang ditunjukkan di bawah ini diambil dari dua distribusi normal yang berbeda, masing-masing dengan mean dan standar deviasi tertentu:

Untuk menghitung perkiraan yang masuk akal dari mean "benar" dan parameter standar deviasi untuk distribusi merah, kita dapat dengan mudah melihat titik merah dan mencatat posisi masing-masing, dan kemudian menggunakan rumus yang sudah dikenal (dan juga untuk kelompok biru) .

Sekarang perhatikan kasus di mana kita tahu bahwa ada dua kelompok titik, tetapi kita tidak bisa melihat titik mana yang termasuk kelompok mana. Dengan kata lain, warnanya disembunyikan:

Sama sekali tidak jelas bagaimana membagi poin menjadi dua kelompok. Kami sekarang tidak dapat hanya melihat posisi dan menghitung estimasi untuk parameter distribusi merah atau distribusi biru.

Di sinilah EM dapat digunakan untuk menyelesaikan masalah.

Menggunakan EM untuk memperkirakan parameter

Berikut adalah kode yang digunakan untuk menghasilkan poin yang ditunjukkan di atas. Anda dapat melihat rata-rata aktual dan standar deviasi dari distribusi normal tempat titik diambil. Variabel reddan bluememegang posisi masing-masing titik dalam kelompok merah dan biru masing-masing:

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

Jika kita dapat melihat warna setiap titik, kita akan mencoba dan memulihkan cara dan standar deviasi menggunakan fungsi pustaka:

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

Tapi karena warnanya tersembunyi dari kami, kami akan memulai proses EM ...

Pertama, kita hanya menebak nilai untuk parameter masing-masing kelompok ( langkah 1 ). Dugaan ini tidak harus bagus:

# 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

Tebakan yang cukup buruk - artinya terlihat jauh dari "tengah" kelompok poin mana pun.

Untuk melanjutkan EM dan meningkatkan tebakan ini, kami menghitung kemungkinan setiap titik data (terlepas dari warna rahasianya) yang muncul di bawah tebakan ini untuk mean dan standar deviasi ( langkah 2 ).

Variabel both_coloursmemegang setiap titik data. Fungsi stats.normmenghitung probabilitas titik di bawah distribusi normal dengan parameter yang diberikan:

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)

Ini memberi tahu kita, misalnya, bahwa dengan perkiraan kami saat ini, titik data di 1,761 jauh lebih mungkin untuk menjadi merah (0,189) daripada biru (0,00003).

Kami dapat mengubah dua nilai kemungkinan ini menjadi bobot ( langkah 3 ) sehingga jumlahnya menjadi 1 sebagai berikut:

likelihood_total = likelihood_of_red + likelihood_of_blue

red_weight = likelihood_of_red / likelihood_total
blue_weight = likelihood_of_blue / likelihood_total

Dengan perkiraan kami saat ini dan bobot kami yang baru dihitung, kami sekarang dapat menghitung estimasi baru, mungkin lebih baik, untuk parameter ( langkah 4 ). Kita membutuhkan fungsi untuk mean dan fungsi untuk deviasi standar:

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)

Ini terlihat sangat mirip dengan fungsi biasa dengan mean dan standar deviasi data. Perbedaannya adalah penggunaan weightparameter yang memberikan bobot pada setiap titik data.

Pembobotan ini adalah kunci EM. Semakin besar bobot warna pada suatu titik data, semakin banyak titik data tersebut mempengaruhi taksiran selanjutnya untuk parameter warna itu. Pada akhirnya, ini memiliki efek menarik setiap parameter ke arah yang benar.

Tebakan baru dihitung dengan fungsi-fungsi ini:

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

Proses EM kemudian diulangi dengan tebakan baru ini dari langkah 2 dan seterusnya. Kita dapat mengulangi langkah-langkah untuk sejumlah iterasi (katakanlah 20), atau sampai kita melihat parameter menyatu.

Setelah lima iterasi, kami melihat dugaan buruk awal kami mulai menjadi lebih baik:

Setelah 20 iterasi, proses EM memiliki lebih atau kurang konvergen:

Sebagai perbandingan, berikut adalah hasil dari proses EM dibandingkan dengan nilai yang dihitung di mana informasi warna tidak disembunyikan:

          | 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

Catatan: jawaban ini diadaptasi dari jawaban saya di Stack Overflow di sini .

Mengikuti jawaban Zhubarb, saya menerapkan contoh “melempar koin” Do and Batzoglou di GNU R. Perhatikan bahwa saya menggunakan mlefungsi stats4paket - ini membantu saya untuk memahami lebih jelas bagaimana EM dan MLE terkait.

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

Semua hal di atas terlihat seperti sumber yang bagus, tetapi saya harus menautkan ke contoh yang bagus ini. Ini menyajikan penjelasan yang sangat sederhana untuk menemukan parameter untuk dua baris dari serangkaian poin. Tutorialnya adalah oleh Yair Weiss saat berada di MIT.

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

Jawaban yang diberikan oleh Zhubarb bagus, tapi sayangnya itu dalam Python. Di bawah ini adalah implementasi Java dari algoritma EM yang dieksekusi pada masalah yang sama (diajukan dalam artikel oleh Do and Batzoglou, 2008). Saya telah menambahkan beberapa printf ke output standar untuk melihat bagaimana parameter bertemu.

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

Kode Java berikut di bawah ini:

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

Baiklah, saya sarankan Anda membaca buku tentang R karya Maria L Rizzo. Salah satu bab berisi penggunaan algoritma EM dengan contoh numerik. Saya ingat membaca kode untuk pemahaman yang lebih baik.

Juga, cobalah untuk melihatnya dari sudut pandang pengelompokan di awal. Kerjakan dengan tangan, masalah pengelompokan di mana 10 pengamatan diambil dari dua kepadatan normal yang berbeda. Ini seharusnya membantu. Ambil bantuan dari 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