### An simulation example
### Author: Lei Pang, Judy Wang & Wenbin Lu
### Department of Statistics, North Carolina State University
### Created Date:  Sep 30, 2010

### Meidan regression for error scenario 1; 15% censoring
### Five methods: IS, BT-boot, S1, S2, S3

rm(list=ls())

source("Functions-IS.R")

N=2000    #number of simulation
n=200   #sample size
tau=0.5 #quantile level
h1=n^-0.1  #bandwidth for kernel-smoothing
h2=n^-0.3
h3=n^-0.7
const1=sqrt(3)/pi #constants for generating logistic and extreme error
const2=sqrt(6)/pi
CP=NULL    #censoring proportion

beta=c(1,1,1) #true parameter

###############  Case1  ##############

set.seed(123456789)
cp=rep(0,N)
Dat=NULL
IS=S1=S2=S3=BTboot=array(0,c(3,4,N))

#generate data
for (i in 1:N)
{
x1 = runif(n,0,4)
x2 = rbinom(n,1,0.5)
x=cbind(1,x1, x2)  #covariates
e = rnorm(n)       #random error;scenario 1
T = as.vector((x%*%beta + e)) #log(t)
k=23.3             #censoring time is unif(0,k)
C = runif(n, 0, k)            #log(c)
y=pmin(T,C)
di=1*(T<=C)                   #delta, censoring indicator
IS[,,i]=IS.func(y, x, di, tau)$Coef  #IS fitting
S1[,,i]=KS.func(y, x, di, tau, h1)$Coef  #kernel-smoothing fitting
S2[,,i]=KS.func(y, x, di, tau, h2)$Coef 
S3[,,i]=KS.func(y, x, di, tau, h3)$Coef 
Dat[[i]]=cbind(y,di,x)        #save data for bootstrap later
cp[i]=1-mean(di)              #record censoring proportions
print(i)
}
mean(cp)                      #show censoring proportion
CP=cbind(CP,c(mean(cp),sd(cp)))

##get seeds for 500 bootstrap resamples
set.seed(37289)
nboot=500
seeds = sample(1:9000000, N)

for (i in 1:N)
{
s=Dat[[i]]
seed=seeds[i]
##get BT-boot estimation
BTboot[,,i]=Boot.func(s[,1], s[,3:5], s[,2], tau, nboot, seed)$Coef[,1:4]
print(i)
}

###process result, generate a table, and write to csv file
table=rbind(processor(IS,beta),processor(BTboot,beta),processor(S1,beta),processor(S1,beta),processor(S3,beta))
colnames(table)=c("Bias","MSE","SE/SD","CovP","Bias","MSE","SE/SD","CovP","Bias","MSE","SE/SD","CovP")
rownames(table)=c("IS","BT-boot","S1","S2","S3")
write.csv(table,"scenario 1;tau=0.5;15%censoring.csv")