library(norm)
library(mvtnorm)
library(car)

sim <- 10000
res <- matrix(0,sim,2)
N <- 500
mean <- c(1.2,0.8)
sigma <- matrix(c(1.1,0.9,0.9,1.4),2,2)
mp <- c(0.7,0.7)

############################################################
############################################################
######################### EM ###############################
############################################################
############################################################


rngseed(1234567)


time_start <- proc.time()[3]

for (i in 1:sim) 
{

mdata <- rmvnorm(N,mean,sigma)
MCAR <- matrix(c(rbinom(N,1,mp[1]),rbinom(N,1,mp[2])),N,2)

mdata[MCAR[,1]==1,1] = NA
mdata[MCAR[,2]==1,2] = NA

s <- prelim.norm(mdata)
mytheta <- em.norm(s,showits=FALSE)

res[i,1] <- getparam.norm(s,mytheta)$mu[1]-getparam.norm(s,mytheta)$mu[2]

if((i/100-floor(i/100))==0) {
print(i);
flush.console();
}

}
proc.time()[3] - time_start

mean(res[!is.nan(res[,1]),1])
sd(res[!is.nan(res[,1]),1])

############################################################
############################################################
######################### Form A       #####################
############################################################
############################################################

rngseed(1234567)

time_start <- proc.time()[3]

for (i in 1:sim) 
{

data <- rmvnorm(N,mean,sigma)
MCAR <- matrix(c(rbinom(N,1,mp[1]),rbinom(N,1,mp[2])),N,2)

data[MCAR[,1]==1,1] = NA
data[MCAR[,2]==1,2] = NA

alldata <- cbind(data,data[,1]-data[,2])

n <- length(alldata[(MCAR[,2]==0)&(MCAR[,1]==0),1])
m1 <- length(alldata[(MCAR[,1]==0)&(MCAR[,2]==1),1])
m2 <- length(alldata[(MCAR[,1]==1)&(MCAR[,2]==0),1])

theta_hat <- mean(alldata[(MCAR[,2]==0)&(MCAR[,1]==0),3])
var_theta_hat <- var(alldata[(MCAR[,2]==0)&(MCAR[,1]==0),3])/(n-1)

beta_hat <- matrix(colMeans(alldata[(MCAR[,1]==0)&(MCAR[,2]==0),1:2]),2,1)
beta_tilde <- matrix(c(mean(alldata[(MCAR[,1]==0)&(MCAR[,2]==1),1]), 
                       mean(alldata[(MCAR[,1]==1)&(MCAR[,2]==0),2])),2,1)

for_inv <- cov(alldata[(MCAR[,2]==0)&(MCAR[,1]==0),1:2])/(n-1)
cov <- cov(alldata[(MCAR[,2]==0)&(MCAR[,1]==0),3],
           alldata[(MCAR[,2]==0)&(MCAR[,1]==0),1:2])/(n-1)
for_inv[1,1] <- for_inv[1,1]*(1+n/m1)
for_inv[2,2] <- for_inv[2,2]*(1+n/m2)

if(det(for_inv)>0) {
inv <- inv(for_inv)
theta <- theta_hat - cov %*% inv %*% (beta_hat - beta_tilde)
var_theta <- var_theta_hat - cov %*% inv %*% t(cov)
}

res[i,1] <- theta
res[i,2] <- theta_hat

if((i/100-floor(i/100))==0) {
print(i);
flush.console();
}

}
proc.time()[3] - time_start

mean(res[!is.nan(res[,1]),1])
sd(res[!is.nan(res[,1]),1])
mean(res[!is.nan(res[,2]),2])
sd(res[!is.nan(res[,2]),2])

############################################################
############################################################
######################### Form B          ##################
############################################################
############################################################

rngseed(1234567)

time_start <- proc.time()[3]

for (i in 1:sim) 
{

data <- rmvnorm(N,mean,sigma)
MCAR <- matrix(c(rbinom(N,1,mp[1]),rbinom(N,1,mp[2])),N,2)

data[MCAR[,1]==1,1] = NA
data[MCAR[,2]==1,2] = NA

alldata <- cbind(data,data[,1]-data[,2])

n       <- length(alldata[(MCAR[,1]==0)&(MCAR[,2]==0),1])
NplusM1 <- length(alldata[(MCAR[,1]==0),1])
NplusM2 <- length(alldata[(MCAR[,2]==0),2])

theta_hat      <- mean(alldata[(MCAR[,1]==0)&(MCAR[,2]==0),3])
var_theta_hat  <-  var(alldata[(MCAR[,1]==0)&(MCAR[,2]==0),3])/(n-1)

beta_hat       <-   matrix(colMeans(alldata[(MCAR[,1]==0)&(MCAR[,2]==0),1:2]),2,1)
beta_tilde_hat <-   matrix(c(mean(alldata[MCAR[,1]==0,1]), 
                             mean(alldata[MCAR[,2]==0,2])),2,1)

var_beta_hat   <- cov(alldata[(MCAR[,1]==0)&(MCAR[,2]==0),1:2])/(n-1)
var_beta_tilde_hat <- c(var(alldata[MCAR[,1]==0,1])/(NplusM1-1), 
                               var(alldata[MCAR[,2]==0,2])/(NplusM2-1))

diag_var_beta_hat <- diag(var_beta_hat)
deltaV <- diag_var_beta_hat-var_beta_tilde_hat
deltaVinv <- deltaV^{-1}

for_inv <- var_beta_hat + 
           diag(diag_var_beta_hat) %*% diag(deltaVinv) %*% diag(diag_var_beta_hat) 
           - diag(diag_var_beta_hat)

cov <- cov(alldata[(MCAR[,2]==0)&(MCAR[,1]==0),3],
           alldata[(MCAR[,2]==0)&(MCAR[,1]==0),1:2])/(n-1)

if(det(for_inv)>0) {
 inv <- inv(for_inv)
 theta <- theta_hat - cov %*% inv %*% diag(diag_var_beta_hat) %*% diag(deltaVinv) %*% (beta_hat - beta_tilde_hat)
 var_theta <- var_theta_hat - cov %*% inv %*% t(cov)
}

res[i,1] <- theta
res[i,2] <- theta_hat

if((i/100-floor(i/100))==0) {
print(i);
flush.console();
}

}
proc.time()[3] - time_start

mean(res[!is.nan(res[,1]),1])
sd(res[!is.nan(res[,1]),1])
mean(res[!is.nan(res[,2]),2])
sd(res[!is.nan(res[,2]),2])