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

############################################################
############################################################
######################### 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

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

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

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

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

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

var_beta_zero1 <- var(boot_res1)
var_beta_zero2 <- var(boot_res2)

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

c1 <- for_inv[1,1] / (for_inv[1,1] - var_beta_zero1)
c2 <- for_inv[2,2] / (for_inv[2,2] - var_beta_zero2)

c <- matrix(c(c1,0,0,c2),2,2);

for_inv[1,1] <- for_inv[1,1]*c1;
for_inv[2,2] <- for_inv[2,2]*c2;

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 %*% (c %*% (beta_hat - beta_zero))
}
}
}
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] 6398.23
# > 
# > mean(res[!is.nan(res[,1]) & res[,1] < 100,1])
# [1] -0.007946366
# > sd(res[!is.nan(res[,1]) & res[,1] < 100,1])
# [1] 0.1508748
# > mean(res[!is.nan(res[,2]) & res[,2] < 100,2])
# [1] -0.4413612
# > sd(res[!is.nan(res[,2]) & res[,2] < 100,2])
# [1] 0.1135784