z <- runif(10000) 
y <- rnorm(10000)+z-2 
ct <- runif(10000)*3-2.5 
failind <- ifelse(ct<y,0,1) 
y <- pmin(y,ct) 
table(failind) 
zbar <- mean(z) 
zc <- z-zbar 
z2 <- sum(zc^2) 
## EM 
b <- c(0,0,1)
bold <- c(-9,-9,-9) 
xb <- b[1]+b[2]*z 
while(sum((bold-b)^2)/sum(b^2) > .00001^2) { 
                                        # E step 
  e0 <- (y-xb)/b[3] 
  t1 <- b[3]*dnorm(e0)/(1-pnorm(e0)) 
  ystar <- failind*y+(1-failind)*(xb+t1) 
  wstar <- failind*y*y+(1-failind)*(b[3]^2+xb^2+t1*(y+xb)) 
                                        # M step 
  bold <- b 
  b[2] <- sum(ystar*zc)/z2 
  b[1] <- mean(ystar)-b[2]*zbar 
  xb <- b[1]+b[2]*z 
  b[3] <- sqrt(sum(wstar-2*ystar*xb+xb^2)/length(y)) 
  print(b) 
 }