#! /bin/Rscript

run_em <- function(data, cutoff = 1.50, theta_init = 0.5, iter = 100, tolerance = 1e-4) {
    missing = data[data == cutoff]
    observed = data[data != cutoff]
    n1 = length(observed)
    n2 = length(missing)
    n = n1 + n2
    theta_current = theta_init
    theta_new = 0
    for (i in 1:iter) {
        theta_current = theta_new
        x_bar = mean(observed)
        theta_new = n1/n * x_bar +
                    n2/n * theta_current +
                    n2/n *
                        dnorm(cutoff - theta_current)/(1 - pnorm(cutoff - theta_current))
        if(theta_new - theta_current <= tolerance) {
            break()
        }
    }
    cat("theta is ", theta_new)
}

data <- c(2.01,0.74,0.68,1.50,1.47,1.50,1.50,1.52,0.07,-0.04,-0.21,0.05,-0.09,0.67,0.14)
run_em(data)