library(norm)
library(car)

sim <- 10000
res <- matrix(0,sim,2)
N <- 2500
mp <- c(0.7,0.7)
alpha <- 1
beta <- 1
gamma <- 0.9
tol <- 0.00000001
boot_size <- 100
imps <- 20;

############################################################
############################################################
######################### median ###########################
############################################################
############################################################

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(0,N,2)
MCAR[,2] <- rbinom(N,1,exp(data[,1])/(1+exp(data[,1])))

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