library(norm)
library(car)
library(Bhat)

sim <- 10000
res <- matrix(0,sim,2)
N <- 500
mp <- c(0.7,0.7)
alpha <- 1
beta <- 1
gamma <- 0.1
tol <- 0.00000001
boot_size <- 100

############################################################
############################################################
######################### MCAR - WLS #######################
############################################################
############################################################

rngseed(1234567)

time_start <- proc.time()[3]

for (i in 1:sim) 
{

data <- matrix(0,nrow=N,ncol=2);
data[,1] <- rexp(n=N,alpha);
data[,2] <- 1-exp(-0.5);

for (j in 1:N) {
tmp <- 1;
u <- runif(1);
while (abs(tmp)>tol) {
F <- 1 - (1+gamma*data[j,2]) * exp(-data[j,2] * (1+gamma*data[j,1]));
f <- ( (1+gamma*data[j,2]) * (1+gamma*data[j,1]) - gamma ) * exp(-data[j,2]*(1+gamma*data[j,1]));
tmp <- (F - u) / f; 
if(is.na(tmp)) { tmp = 0; j = j - 1;}
data[j,2] <- data[j,2] - tmp;
}
}

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

theta_hat <- median(alldata[(MCAR[,1]==0)&(MCAR[,2]==0),3])

beta_hat <- matrix(0,2,1)
beta_hat[1] <- median(alldata[(MCAR[,1]==0)&(MCAR[,2]==0),1])
beta_hat[2] <- median(alldata[(MCAR[,1]==0)&(MCAR[,2]==0),2])

beta_tilde <- matrix(0,2,1)
beta_tilde[1] <- median(alldata[(MCAR[,1]==0)&(MCAR[,2]==1),1])
beta_tilde[2] <- median(alldata[(MCAR[,1]==1)&(MCAR[,2]==0),2])

for_boot <- alldata[(MCAR[,1]==0)&(MCAR[,2]==0),1:2];
boot_res <- matrix(0,boot_size,2);
for(boot in 1:boot_size) {
tmp <- for_boot[sample(1:n,replace=TRUE),1:2];
boot_res[boot,1] <- median(tmp[,1])
boot_res[boot,2] <- median(tmp[,2])
}

for_inv <- cov(boot_res)

cov <- matrix(c(for_inv[1,1] - for_inv[1,2], for_inv[1,2] - for_inv[2,2]),1,2);

for_inv[1,1] <- for_inv[1,1]*(1+n/m1)
for_inv[2,2] <- for_inv[2,2]*(1+n/m2)

res[i,1] <- theta_hat;
if (!is.na(max(for_inv))) {
if(max(for_inv)<10000000) {
if(det(for_inv)>tol) {
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]) & res[,1] < 100,1])
sd(res[!is.nan(res[,1]) & res[,1] < 100,1])
mean(res[!is.nan(res[,2]) & res[,2] < 100,2])
sd(res[!is.nan(res[,2]) & res[,2] < 100,2])

############################################################
############################################################
######################### MCAR - MI ########################
############################################################
############################################################

rngseed(1234567)

time_start <- proc.time()[3]

for (i in 1:sim) 
{

data <- matrix(0,nrow=N,ncol=2);
data[,1] <- rexp(n=N,alpha);
data[,2] <- 1-exp(-0.5);

for (j in 1:N) {
tmp <- 1;
u <- runif(1);
while (abs(tmp)>tol) {
F <- 1 - (1+gamma*data[j,2]) * exp(-data[j,2] * (1+gamma*data[j,1]));
f <- ( (1+gamma*data[j,2]) * (1+gamma*data[j,1]) - gamma ) * exp(-data[j,2]*(1+gamma*data[j,1]));
tmp <- (F - u) / f; 
if(is.na(tmp)) { tmp = 0; j = j - 1;}
data[j,2] <- data[j,2] - tmp;
}
}

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

s <- prelim.norm(data)
thetahat <- em.norm(s,showits=FALSE)

med <- 0;
for (t in 1:imps) {
 ximp <- imp.norm(s,thetahat,data);
 med <- med + median(ximp[,1]) - median(ximp[,2])
}
res[i,1] <- med / imps;

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