"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) } } "em.h"<- function(x) { dnorm(x)/(1 - pnorm(x)) }