"em.mc"<-
function(x, m, th.init)
{
        th <- th.init
        for(i in 1:length(m)) {
                p <- th/(2 + th)
                z <- rbinom(m[i], 125, p)
                me <- mean(z)
                th <- (me + x[4])/(me + x[2] + x[3] + x[4])
                cat(i, round(th, 5), fill = T)
        }
}
"em.reg"<-
function(x, n)
{
 m_dim(x)[1]
        a <- lm(x[, 2] ~ x[, 1])
        b <- coef(a)
        sig <- sqrt(deviance(a)/(m-2))
        cat(b, sig, fill = T)
        for(i in 1:n) {
                mu <- b[1] + b[2] * x[, 1]
                ez <- mu + sig * em.h((x[, 2] - mu)/sig)
                ez[x[,3]==0] <- x[x[,3]==0, 2]
                a <- lm(ez ~ x[, 1])
                b <- coef(a)
                ss <- sum((ez[x[,3]==0] - mu[x[,3]==0])^2)/m
                ss <- ss + (sig^2 * sum(((x[x[,3]==1, 2] - mu[x[,3]==1])/sig) * em.h((
                        x[x[,3]==1, 2] - mu[x[,3]==1])/sig) + 1))/m
                sig <- sqrt(ss)
                cat(b, sig, fill = T)
                abline(a, col = i)
        }
}