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.5
tol <- 0.00000001
boot_size <- 100

############################################################
############################################################
######################### means ############################
############################################################
############################################################

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

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

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

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 <- for_inv[1,1] - for_inv[1,2];
for_inv[1,1] <- for_inv[1,1]*(1+n/m1)

res[i,1] <- theta_hat;
for_inv1 <- for_inv[1,1];

if (!is.na(for_inv1)) {
if(for_inv1<10000000) {
if(for_inv1>tol) {
theta <- theta_hat - cov * (beta_hat - beta_tilde) / for_inv1
}
}
}
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])

# > proc.time()[3] - time_start
# [1] 5139.93
# > 
# > mean(res[!is.nan(res[,1]) & res[,1] < 100,1])
# [1] 0.0622531
# > sd(res[!is.nan(res[,1]) & res[,1] < 100,1])
# [1] 0.1508747
# > mean(res[!is.nan(res[,2]) & res[,2] < 100,2])
# [1] -0.4399997
# > sd(res[!is.nan(res[,2]) & res[,2] < 100,2])
# [1] 0.1144867
# >