## (X, Y) ~ BVN(mu_x, mu_y, sigma_xx, sigma_yy, sigma_xy).  Some X, Y missing

## True values: means = 0, variances = 1, correlation = 0.5

Sigma <- matrix(c(1, 0.5, 0.5, 1), 2, 2)
Tmat <- t(chol(Sigma))

xy <- Tmat %*% matrix(rnorm(2000), nrow = 2)
xy <- as.data.frame(t(xy))
colnames(xy) <- c("x", "y")

## Set some values to NA

xyna <- xy
f <- 0.3
i <- sample(seq_len(nrow(xyna)), round(f * nrow(xyna)))
w <- ifelse(runif(length(i)) < 0.5, "x", "y")
xyna$x[i[w == "x"]] <- NA
xyna$y[i[w == "y"]] <- NA



naiveBVN <- function(x, y, niter = 10)
{
    nax <- is.na(x)
    nay <- is.na(y)
    nona <- !(nax | nay)
    mux <- mean(x[!nax])
    muy <- mean(y[!nay])
    sigma <- var(cbind(x[nona], y[nona]))
    for (i in seq_len(niter))
    {
        print(unname(diag(sigma)))
        ## impute based on current estimate
        xx <- x
        yy <- y
        xx[nax] <- mux + sigma[1, 2] / sigma[2, 2] * (yy[nax] - muy)
        yy[nay] <- muy + sigma[1, 2] / sigma[1, 1] * (xx[nay] - mux)
        ## plot(xx, yy, pch = 16)
        mux <- mean(xx)
        muy <- mean(yy)
        sigma <- var(cbind(xx, yy))
    }
    list(mux = mux, muy = muy, sigma = sigma)
}

naiveBVN(xyna$x, xyna$y, niter = 20)