2- You may have question marks in your head, especially
regarding where the probabilities in the Expectation step come from. Please have
a look at the explanations on this maths stack exchange page.

3- Look at/run this code that I wrote in Python that
simulates the solution to the coin-toss problem in the EM tutorial paper of item
1:

P.S The code may be suboptimal, but it does the job.

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

Expectation-Maximization in CSharp

Jump to: navigation,
search

This example requires Emgu
CV 1.5.0.0

using System.Drawing;

using Emgu.CV.Structure;

using Emgu.CV.ML;

using Emgu.CV.ML.Structure;

...
int N = 4; //number of clusters

int N1 = (int)Math.Sqrt((double)4);
Bgr[] colors = new Bgr[] {

   new Bgr(0, 0, 255),

   new Bgr(0, 255, 0),

   new Bgr(0, 255, 255),

   new Bgr(255, 255, 0)};
int nSamples = 100;
Matrix<float> samples = new Matrix<float>(nSamples, 2);

Matrix<Int32> labels = new Matrix<int>(nSamples, 1);

Image<Bgr, Byte> img = new Image<Bgr,byte>(500, 500);

Matrix<float> sample = new Matrix<float>(1, 2);
CvInvoke.cvReshape(samples.Ptr, samples.Ptr, 2, 0);

for (int i = 0; i < N; i++)

{

   Matrix<float> rows = samples.GetRows(i * nSamples / N, (i + 1) * nSamples / N, 1);

   double scale = ((i % N1) + 1.0) / (N1 + 1);

   MCvScalar mean = new MCvScalar(scale * img.Width, scale * img.Height);

   MCvScalar sigma = new MCvScalar(30, 30);

   ulong seed = (ulong)DateTime.Now.Ticks;

   CvInvoke.cvRandArr(ref seed, rows.Ptr, Emgu.CV.CvEnum.RAND_TYPE.CV_RAND_NORMAL, mean, sigma);

}

CvInvoke.cvReshape(samples.Ptr, samples.Ptr, 1, 0);
using (EM emModel1 = new EM())

using (EM emModel2 = new EM())

{

   EMParams parameters1 = new EMParams();

   parameters1.Nclusters = N;

   parameters1.CovMatType = Emgu.CV.ML.MlEnum.EM_COVARIAN_MATRIX_TYPE.COV_MAT_DIAGONAL;

   parameters1.StartStep = Emgu.CV.ML.MlEnum.EM_INIT_STEP_TYPE.START_AUTO_STEP;

   parameters1.TermCrit = new MCvTermCriteria(10, 0.01);

   emModel1.Train(samples, null, parameters1, labels);
EMParams parameters2 = new EMParams();

   parameters2.Nclusters = N;

   parameters2.CovMatType = Emgu.CV.ML.MlEnum.EM_COVARIAN_MATRIX_TYPE.COV_MAT_GENERIC;

   parameters2.StartStep = Emgu.CV.ML.MlEnum.EM_INIT_STEP_TYPE.START_E_STEP;

   parameters2.TermCrit = new MCvTermCriteria(100, 1.0e-6);

   parameters2.Means = emModel1.GetMeans();

   parameters2.Covs = emModel1.GetCovariances();

   parameters2.Weights = emModel1.GetWeights();
emModel2.Train(samples, null, parameters2, labels);
#region Classify every image pixel

   for (int i = 0; i < img.Height; i++)

      for (int j = 0; j < img.Width; j++)

      {

         sample.Data[0, 0] = i;

         sample.Data[0, 1] = j;

         int response = (int) emModel2.Predict(sample, null);
Bgr color = colors[response];
img.Draw(

            new CircleF(new PointF(i, j), 1),

            new Bgr(color.Blue*0.5, color.Green * 0.5, color.Red * 0.5 ),

            0);

      }

   #endregion
#region draw the clustered samples

   for (int i = 0; i < nSamples; i++)

   {

      img.Draw(new CircleF(new PointF(samples.Data[i, 0], samples.Data[i, 1]), 1), colors[labels.Data[i, 0]], 0);

   }

   #endregion
Emgu.CV.UI.ImageViewer.Show(img);

}

ExpectationMaximum