# 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)
n <- length(Time);
par(mfrow=c(3,3))
   
# premier graphe : distribution de la variable Time
hist(Time, nclass=100, xlab='Population Distribution');

# deuxième graphe : on estime la moyenne sur un échantillon de k=20
# on estime par Monte-Carlo la (vraie) loi de cet estimateur
# et on montre un IC empirique
N <- 10000;
k <-20; # taille de l'échantillon sur lequel on travaille

Ave<-array(0, N);
for (i in 1:N){
  x<-sample(Time, k, replace=F);
  Ave[i]<-mean(x);
}
hist(Ave, nclass=100, xlim=c(0,30), xlab='Sample Distribution');

mu <- mean(Time); sAve<- sort(Ave)
alpha=5/100; 
IC_Ave = c(sAve[floor(alpha/2*N)], sAve[ceiling((1-alpha/2)*N)]);
lines(IC_Ave, c(-1,-1), col='green', lwd=3)
points(mu,0, col='red', pch='x', cex=2)

# troisième graphe : QQ plot pour voir si la loi de l'estimateur est gaussienne
qqnorm(Ave, xlab='Is it Gaussian?')

# premier graphe bis : on fixe un échantillon de k=20 
# => sa loi empirique ressemble à peu près à la loi de Time
s1 <- sample(Time, k, replace=F);
hist(s1, nclass=k, xlab='Sample 1');


# deuxième graphe bis, et ter : on estime la moyenne sur l'échantillon
# fixé de de k=20 et on montre les IC bootstrap
BootAve<-array(0, N);
for (i in 1:N){
  x<-sample(s1, replace=T);
  bootAve[i]<-mean(x);
}
hist(bootAve, nclass=100, xlim=c(0,30), xlab='Sample 1, Bootstrap 1');
xb <- mean(bootAve); sbAve<- sort(bootAve)
alpha=5/100; 
IC_std = mean(s1) + qt(1-alpha/2, k-1)*sd(s1)/sqrt(k)*c(-1,1);
IC_Ave = c(sbAve[floor(alpha/2*N)], sbAve[ceiling((1-alpha/2)*N)]);
lines(IC_Ave, c(2,2), col='blue', lwd=3)
lines(IC_std, c(-4,-4), col='green', lwd=3)
points(xb,0, col='red', pch='x', cex=2)


bootAve<-array(0, N);
for (i in 1:N){
  x<-sample(s1, replace=T);
  bootAve[i]<-mean(x);
}
hist(bootAve, nclass=100, xlim=c(0,30), xlab='Sample 1, Bootstrap 2');
xb <- mean(bootAve); sbAve<- sort(bootAve)
alpha=5/100; 
IC_std = mean(s1) + qt(1-alpha/2, k-1)*sd(s1)/sqrt(k)*c(-1,1);
IC_Ave = c(sbAve[floor(alpha/2*N)], sbAve[ceiling((1-alpha/2)*N)]);
lines(IC_Ave, c(2,2), col='blue', lwd=3)
lines(IC_std, c(-4,-4), col='green', lwd=3)
points(xb,0, col='red', pch='x', cex=2)

# on recommence avec un nouvel échantillon de k=20

s2 <- sample(Time, k, replace=F);
hist(s2, nclass=k, xlab='Sample 2');

bootAve<-array(0, N);
for (i in 1:N){
  x<-sample(s2, replace=T);
  bootAve[i]<-mean(x);
}
hist(bootAve, nclass=100, xlim=c(0,30), xlab='Sample 2, Bootstrap 1');
xb <- mean(bootAve); sbAve<- sort(bootAve)
alpha=5/100; 
IC_std = mean(s2) + qt(1-alpha/2, k-1)*sd(s2)/sqrt(k)*c(-1,1);
IC_Ave = c(sbAve[floor(alpha/2*N)], sbAve[ceiling((1-alpha/2)*N)]);
lines(IC_Ave, c(2,2), col='blue', lwd=3)
lines(IC_std, c(-4,-4), col='green', lwd=3)
points(xb,0, col='red', pch='x', cex=2)

bootAve<-array(0, N);
for (i in 1:N){
  x<-sample(s2, replace=T);
  bootAve[i]<-mean(x);
}
hist(bootAve, nclass=100, xlim=c(0,30), xlab='Sample 2, Bootstrap 2');
xb <- mean(bootAve); sbAve<- sort(bootAve)
alpha=5/100; 
IC_std = mean(s2) + qt(1-alpha/2, k-1)*sd(s2)/sqrt(k)*c(-1,1);
IC_Ave = c(sbAve[floor(alpha/2*N)], sbAve[ceiling((1-alpha/2)*N)]);
lines(IC_Ave, c(2,2), col='blue', lwd=3)
lines(IC_std, c(-4,-4), col='green', lwd=3)
points(xb,0, col='red', pch='x', cex=2)