Basic Expectation Maximization example

### assign probability of getting head to two unfair coins

(head.prob1.init <- runif(1))
(head.prob2.init <- runif(1))

### simulate coin tosses

k <- 10   # coin selections; article: 5
n <- 100  # coin tosses per selection; article: 10

coin  <- numeric(k) # coin identify
heads <- numeric(k) # numer of heads in a turn

for (i in 1:k) {

   coin[i]   <- ifelse(runif(1)<=0.5,1,2)

   tosses    <- rbinom(n, 1, ifelse(coin[i]==1,head.prob1.init,head.prob2.init))

   heads[i]  <- length(which(tosses==1))

}

### esitmate coins probabilities with MLE

head.prob1.est <- sum(heads[which(coin==1)]) / (length(which(coin==1))*n)
head.prob2.est <- sum(heads[which(coin==2)]) / (length(which(coin==2))*n)

# comparison of MLE estimates to the actual probabilties

head.prob1.init; head.prob1.est
head.prob2.init; head.prob2.est

### esitmate coins probabilities with MaximumExpectation

heads2 <- heads # used in the example used in the article: c(5,9,8,4,7)

prob1.maxexp <- runif(1) # article: 0.6
prob2.maxexp <- runif(1) # article: 0.5

prev.prob1 <- 0
prev.prob2 <- 0

cp <- 3 # rounding precision

while ( round(prob1.maxexp,cp)!=round(prev.prob1,cp) &&
        round(prob2.maxexp,cp)!=round(prev.prob2,cp) ) {

  prev.prob1 <- prob1.maxexp
  prev.prob2 <- prob2.maxexp

  den1 <- dbinom(heads,n,prob1.maxexp)
  den2 <- dbinom(heads,n,prob2.maxexp)

  h1 <- den1/(den1+den2)*heads
  h2 <- den2/(den1+den2)*heads

  t1 <- den1/(den1+den2)*(n-heads)
  t2 <- den2/(den1+den2)*(n-heads)

  (prob1.maxexp <- sum(h1)/sum(c(h1,t1)))
  (prob2.maxexp <- sum(h2)/sum(c(h2,t2)))

}

heads
coin

# actual probabilities
head.prob1.init
head.prob2.init

# MLE estimates
head.prob1.est
head.prob2.est

# Maximum Expectation estimates
prob1.maxexp
prob2.maxexp