# UE7 Cours VII.
# Utilisation de la fonction boot (de la librairie boot)
# ----------------------------------------------------------------------

library(boot)

# Un premier exemple (fictif) dans lequel on calcule simplement 
# une moyenne et son erreur standard

donnees <- c(0,10,15,8,9,14.2,20,7,3,12,10,6,11.5)
mean(donnees)

# calcul de l'erreur standard (SD)  à partir des données (estimation usuelle)
sqrt(var(donnees))/sqrt(length(donnees))


#maintenant voyons cela par bootstrap

# la fonction à calculer; les indices servent à échantillonner
resultat=function(valeurs,indices){
mean(valeurs[indices])
}

b<-boot(data = donnees, resultat, R = 100, stype = "i")

#les résultats des calculs
ls(b)                 # que contient b ?
#print(c(b$t))         # le vecteur des 1000 ratios calculés
cat( 'Comparaison des moyennes donnes <=> bootstrap : ', mean(donnees),' <=> ', mean(b$t),"\t")
print(b)
plot(b)
p <- c(0.025, 0.975)
quantile(c(b$t),p) 
  
#------------------------------------------------------------------------#
# les donnees sont les population en 1920 et 1930 des 10 villes
# les plus importantes des USA
# on veut calculer le rapport des moyennes des populations
# pour trouver une sorte de coefficient multiplicateur de la
# population de 1920
#------------------------------------------------------------------------#

data(city))
cat( "ratio moyen : ",(r = sum(city$x )/sum(city$u )))

# mais cela ne dit rien sur la précision de cette valeur !
# Pour calculer cette précision on applique la technique du bootstrapping

# la fonction à calculer.  w indique qu'il s'agit de comptage (weight)
ratio <- function(d, w) sum(d$x * w)/sum(d$u * w)

# le ré-échantllonnage

b<-boot(data = city, ratio, R = 999, stype = "w")

# les résultats
print(b)
print(c(b$t))          # le vecteur des 1000 ratios calculés

# la distribution semble suivre une loi exponentielle et non normale :
plot(b)                # un qqnorm plot: qqnorm(b$t); qqline(b$t)

#quelques quantiles
quantile(c(b$t))

# les bornes de l'intervalle de confiance est donné par
p <- c(0.025, 0.975)
quantile(c(b$t),p)