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.2,0.2,1.4),2,2)
mp <- c(0,0.7)

############################################################
############################################################
######################### 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[,2]==1,2] = NA

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

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

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

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

inv <- var(alldata[(MCAR[,1]==0)&(MCAR[,2]==0),1])
cov <- cov(alldata[(MCAR[,1]==0)&(MCAR[,2]==0),3], alldata[(MCAR[,1]==0)&(MCAR[,2]==0),1])

if(inv > 0) {
theta <- theta_hat - m1 * cov * (beta_hat - beta_tilde) / (inv * (n+m1) )
}

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[,2]==1,2] = NA

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

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

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 <- mean(alldata[(MCAR[,1]==0)&(MCAR[,2]==0),1])
beta_tilde_hat <- mean(alldata[(MCAR[,1]==0),1])

inv <- var(alldata[(MCAR[,1]==0)&(MCAR[,2]==0),1])
cov <- cov(alldata[(MCAR[,1]==0)&(MCAR[,2]==0),3], 
           alldata[(MCAR[,1]==0)&(MCAR[,2]==0),1])

if(inv > 0) {
theta <- theta_hat - cov * (beta_hat - beta_tilde_hat) / inv
}

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])



############################################################
######################### 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[,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])