• API
• FAQ
• Tools
• Archive
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.
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
14.
15. for (i in 1:k) {
16.
17.    coin[i]   <- ifelse(runif(1)<=0.5,1,2)
18.
20.
22.
23. }
24.
25. ### esitmate coins probabilities with MLE
26.
29.
30. # comparison of MLE estimates to the actual probabilties
31.
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.
55.
58.
61.
62.   (prob1.maxexp <- sum(h1)/sum(c(h1,t1)))
63.   (prob2.maxexp <- sum(h2)/sum(c(h2,t2)))
64.
65. }
66.
68. coin
69.
70. # actual probabilities