how about this: http://en.wikibooks.org/wiki/Data_Mining_Algorithms_In_R/Clustering/Expectation_Maximization_(EM)#A_simple_example
I also wrote a simple example in (edit) R a year ago, unfortunately I cannot find it. I will try again to find it later.
EDIT: Here it is -
EM <- function ()
{
### Read file, get necessary cols
dataFile <- read.csv ("wine.csv", head = FALSE, sep = ",")
sl <- dataFile [, 2]
#sw <- dataFile [, 3]
#pl <- dataFile [, 3]
#pw <- dataFile [, 4]
class <- dataFile [, 5]
N <- length (sl)
pi1 <- 0.5
### Init ###
rand1 <- floor (runif (1) * N)
rand2 <- floor (runif (1) * N)
mu1 <- sl [rand1]
mu2 <- sl [rand2]
mean1 <- sum (sl) / N
sigma1 <- sum ((sl - mean1) ** 2) / N
sigma2 <- sigma1
print (mu1)
print (mu2)
print (sigma1)
print (sigma2)
COUNTLIM <- 10
count <- 1
prevmu1 <- 0.0;
prevmu2 <- 0.0;
prevsigma1 <- 0.0;
prevsigma2 <- 0.0;
gamma <- array (0, length (sl))
while (count <= COUNTLIM)
{
gamma <- pi1 * dnorm (sl, mu2, sigma2) / ((1 - pi1) * dnorm (sl, mu1, sigma1) + pi1 * dnorm (sl, mu2, sigma2))
mu1 <- sum ((1 - gamma) * sl) / sum (1 - gamma)
mu2 <- sum ((gamma) * sl) / sum (gamma)
sigma1 <- sum ((1 - gamma) * (sl - mu1) ** 2) / sum (1 - gamma)
sigma2 <- sum ((gamma) * (sl - mu2) ** 2) / sum (gamma)
pi1 <- sum (gamma) / N
print (c (mu1, mu2, sigma1, sigma2, pi1))
if (count == 1)
{
prevmu1 <- mu1;
prevmu2 <- mu2;
prevsigma1 <- sigma1;
prevsigma2 <- sigma2;
}
else
{
val <- ((prevmu1 - mu1) * 2 + (prevmu2 - mu2) * 2 + (prevsigma1 - sigma1) * 2 + (prevsigma2 - sigma2) * 2) ** 0.5;
print (c ("val:", val))
if (val <= 1)
{
break;
}
}
count <- count + 1
}
print (mu1)
print (mu2)
print (sigma1)
print (sigma2)
}
source share