Esempio numerico dell'algoritmo di massimizzazione delle aspettative

Question

Qualcuno potrebbe fornire un semplice esempio numerico dell'algoritmo EM in quanto non sono sicuro delle formule fornite? Uno veramente semplice con 4 o 5 coordinate cartesiane farebbe perfettamente.

Answer 1

che dire di questo: http://en.wikibooks.org/wiki/Data_Mining_Algorithms_In_R/Clustering/Expectation_Maximization_(EM)#A_simple_example

avevo anche scritto un semplice esempio in (modifica) R un anno fa, purtroppo non sono in grado di individuarlo. Proverò di nuovo a trovarlo più tardi.

EDIT: qui è -

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