EM Algorithm

> data(faithful)
> attach(faithful)
        The following object(s) are masked from faithful ( position 3 ) :

         eruptions waiting

    The following object(s) are masked from faithful ( position 4 ) :


         eruptions waiting

    The following object(s) are masked from faithful ( position 5 ) :


         eruptions waiting
> W <- waiting
> hist(W, breaks = 100)
> s <- c(0.5, 40, 90, 16, 16)
> em <- function(W, s) {
+     Ep <- s[1] * dnorm(W, s[2], sqrt(s[4]))/(s[1] * dnorm(W,
+         s[2], sqrt(s[4])) + (1 - s[1]) * dnorm(W, s[3], sqrt(s[5])))
+     s[1] <- mean(Ep)
+     s[2] <- sum(Ep * W)/sum(Ep)
+     s[3] <- sum((1 - Ep) * W)/sum(1 - Ep)
+     s[4] <- sum(Ep * (W - s[2])^2)/sum(Ep)
+     s[5] <- sum((1 - Ep) * (W - s[3])^2)/sum(1 - Ep)
+     return(s)
+ }
> iter <- function(W, s) {
+     s1 <- em(W, s)
+     for (i in 1:5) {
+         if (abs(s[i] - s1[i]) > 1e-04) {
+             s <- s1
+             iter(W, s)
+         }
+         else {
+             s1
+         }
+     }
+     return(s1)
+ }
> p <- iter(W, s)
> p1 <- p[1]
> p2 <- p[2]
> p3 <- p[3]
> p4 <- p[4]
> p5 <- p[5]
> Boot <- function(B) {
+     r <- 0
+     k <- 0
+     bootmean1 <- rep(0, B)
+     bootvar1 <- rep(0, B)
+     bootmean2 <- rep(0, B)
+     bootvar2 <- rep(0, B)
+     for (i in 1:B) {
+         p <- runif(1, 0, 1)
+         if (p < p1) {
+             boot1 <- rnorm(p1 * 272, p2, sqrt(p4))
+             bootmean1[i] <- mean(boot1)
+             bootvar1[i] <- var(boot1)
+             r <- r + 1
+         }
+         else {
+             boot2 <- rnorm((1 - p1) * 272, p3, sqrt(p5))
+             bootmean2[i] <- mean(boot2)
+             bootvar2[i] <- var(boot2)
+             k <- k + 1
+         }
+     }
+     meanbootm1 <- sum(bootmean1)/r
+     meanbootvar1 <- sum(bootvar1)/r
+     meanbootm2 <- sum(bootmean2)/k
+     meanbootvar2 <- sum(bootvar2)/k
+     list(meanbootm1 = meanbootm1, meanbootvar1 = meanbootvar1,
+         meanbootm2 = meanbootm2, meanbootvar2 = meanbootvar2)
+ }
> Boot(1000)
$meanbootm1
[1] 54.232

$meanbootvar1
[1] 29.98183

$meanbootm2
[1] 79.9182

$meanbootvar2
[1] 36.00942