## EM for simple regression
 
## data
 
x<-rnorm(100)
y<-x+.5*rnorm(100)
 
summary(lm(y~x-1))
 
y[91:100]<-NA
 
a0<-1
s0<-.5
 
e<-1
 
while (e>.00001){
 
sx<-sum(x^2)
sy<-sum(y[1:90]^2)+a0^2*sum(x[91:100]^2)+10*s0
sxy<-sum(x[1:90]*y[1:90])+a0*sum(x[91:100]^2)
 
a<-sxy/sx
s<-(sx*sy-sxy^2)/(sx^2)
 
e<-(a0-a)^2+(s0-s)^2
 
a0<-a
s0<-s
}