SHARE
TWEET

Basic Expectation Maximization example

mjaniec Sep 26th, 2013 761 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1.  
  2. ### assign probability of getting head to two unfair coins
  3.  
  4. (head.prob1.init <- runif(1))
  5. (head.prob2.init <- runif(1))
  6.  
  7. ### simulate coin tosses
  8.  
  9. k <- 10   # coin selections; article: 5
  10. n <- 100  # coin tosses per selection; article: 10
  11.  
  12. coin  <- numeric(k) # coin identify
  13. heads <- numeric(k) # numer of heads in a turn
  14.  
  15. for (i in 1:k) {
  16.  
  17.    coin[i]   <- ifelse(runif(1)<=0.5,1,2)
  18.  
  19.    tosses    <- rbinom(n, 1, ifelse(coin[i]==1,head.prob1.init,head.prob2.init))
  20.    
  21.    heads[i]  <- length(which(tosses==1))
  22.  
  23. }
  24.  
  25. ### esitmate coins probabilities with MLE
  26.  
  27. head.prob1.est <- sum(heads[which(coin==1)]) / (length(which(coin==1))*n)
  28. head.prob2.est <- sum(heads[which(coin==2)]) / (length(which(coin==2))*n)
  29.  
  30. # comparison of MLE estimates to the actual probabilties
  31.  
  32. head.prob1.init; head.prob1.est
  33. head.prob2.init; head.prob2.est
  34.  
  35. ### esitmate coins probabilities with MaximumExpectation
  36.  
  37. heads2 <- heads # used in the example used in the article: c(5,9,8,4,7)
  38.  
  39. prob1.maxexp <- runif(1) # article: 0.6
  40. prob2.maxexp <- runif(1) # article: 0.5
  41.  
  42. prev.prob1 <- 0
  43. prev.prob2 <- 0
  44.  
  45. cp <- 3 # rounding precision
  46.  
  47. while ( round(prob1.maxexp,cp)!=round(prev.prob1,cp) &&
  48.         round(prob2.maxexp,cp)!=round(prev.prob2,cp) ) {
  49.  
  50.   prev.prob1 <- prob1.maxexp
  51.   prev.prob2 <- prob2.maxexp
  52.  
  53.   den1 <- dbinom(heads,n,prob1.maxexp)
  54.   den2 <- dbinom(heads,n,prob2.maxexp)
  55.  
  56.   h1 <- den1/(den1+den2)*heads
  57.   h2 <- den2/(den1+den2)*heads
  58.  
  59.   t1 <- den1/(den1+den2)*(n-heads)
  60.   t2 <- den2/(den1+den2)*(n-heads)
  61.  
  62.   (prob1.maxexp <- sum(h1)/sum(c(h1,t1)))
  63.   (prob2.maxexp <- sum(h2)/sum(c(h2,t2)))
  64.  
  65. }
  66.  
  67. heads
  68. coin
  69.  
  70. # actual probabilities
  71. head.prob1.init
  72. head.prob2.init
  73.  
  74. # MLE estimates
  75. head.prob1.est
  76. head.prob2.est
  77.  
  78. # Maximum Expectation estimates
  79. prob1.maxexp
  80. prob2.maxexp
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top