# on efface tout et on charge les données verizon
# variable Time : temps d'attente avant traitement de la demande
# variable Group : type de client (clients directs, clients couverts par un accord)
rm(list=ls())
verizon <- read.table("verizon.txt", header=T, sep="\t");
attach(verizon); names(verizon)

# on affiche un histogramme des données pour se convaincre qu'elles ne sont 
# pas du tout gaussiennes 
par(mfrow=c(2,1));
hist(Time, nclass=100);
qqnorm(Time)

# on bootstrap la moyenne de Time pour vérifier que son estimée, par contre, 
# est bien gaussienne
n <- length(Time);

N<-10000;
bootAve<-array(0, N);
for (i in 1:N){
  x<-sample(Time, replace=T);
  bootAve[i]<-mean(x);
}
par(mfrow=c(2,1))
hist(bootAve, nclass=100);
#plot(sort(mV), (1:N)/N, type="s");
mVs <- sort(bootAve);
alpha=5/100;
IC_std = mean(bootAve) + qt(1-alpha/2, n-1)*sd(Time)/sqrt(n)*c(-1,1); IC_std
IC_boot = c(mVs[floor(alpha/2*N)], mVs[ceiling((1-alpha/2)*N)]); IC_boot
lines(IC_std, c(0,0), col='red', lwd=3)
lines(IC_boot, c(-5,-5), col='green', lwd=3)
points(c(mean(Time), mean(bootAve)), c(0,0), col='black', pch='x', cex=2)
qqnorm(bootAve)

# on définit un fonction d'estimation de densité pour la suite...
epanechnikov <- function (t) {3/4*(1-t^2)*(t>-1)*(t<1)}
plotKernel = function(y, lambda,kernelType, ...){
  k<-length(y);
  tt <- seq(-5, 180, 0.01);
  fh <- tt-tt; 
  for (j in 1:k){
      if (kernelType=='Epanechnikov'){
      contrib <- 1/(k*lambda) * epanechnikov((tt-y[j])/lambda);
    }
      else{ # Gaussian Kernel
      contrib <- 1/(k*lambda) * dnorm((tt-y[j])/lambda);
    }
    fh <- fh + contrib;
    }
  lines(tt, fh, type='l', ...);
  }

# on représente un estimé non paramétrique de la distribution Time pour chaque
# type de clients
par(mfrow=c(1,1))
#hist(Time[Group=='ILEC'], nclass=100)#, col="black")
#hist(Time[Group=='CLEC'], nclass=100)#, col="red")
plot(0,0,type='n', xlim=c(min(Time), max(Time)), ylim=c(0,0.2), xlab='red: clients, blue: agreement', ylab='density');
y<- Time[Group=='ILEC'];
plotKernel(y, 0.5, 'Gaussian', lwd=3, col='red');#'Epanechnikov');
y<- Time[Group=='CLEC'];
plotKernel(y, 2, 'Gaussian', lwd=3, col='blue');

# on fait un test bootstrap pour la différence de temps d'attente moyen
# on dessine les intervalles de confiance de Student et de bootstrap pour la 
# différence entre les deux moyennes
N<-10000;
bootDiff<-array(0, N);
for (i in 1:N){
  x<-sample(Time[Group=='ILEC'], replace=T);
  y<-sample(Time[Group=='CLEC'], replace=T);
  bootDiff[i]<-mean(x)-mean(y);
}
par(mfrow=c(1,1))
hist(bootDiff, nclass=100);
mVs <- sort(bootDiff);
alpha=5/100;
IC_boot = c(mVs[floor(alpha/2*N)], mVs[ceiling((1-alpha/2)*N)]); IC_boot
lines(IC_boot, c(-1,-1), col='red', lwd=3)
ttest <- t.test(Time[Group=='ILEC'], Time[Group=='CLEC'])
lines(ttest$conf.int, c(-5,-5), col='green', lwd=3)
t<-seq(-30, 10, 0.01); n<-length(Time)
m1 <- mean(Time[Group=='ILEC']); n1<- length(Time[Group=='ILEC']);
m2 <- mean(Time[Group=='CLEC']); n2 <- length(Time[Group=='CLEC'])
s <- sqrt((sum((Time[Group=='ILEC']-m1)^2) + sum((Time[Group=='CLEC']-m2)^2))/(n-2))
#--
lines(t, n*dt(sqrt(n1*n2/(n1+n2))*(t-(m1-m2))/s, n-2), col='cyan', lwd=2)