############### EXERCICE 1 

###########################"
#   UN PREMIER EXEMPLE POUR CALCULER LES ESTIMATIONS BOOSTRAP
#  POUR LES COEFFICIENTS d'UNE REGRESSION LINEAIRE
##########################################" 
library(boot)

data(dogs)
#attach(dogs)
m1.lm=lm(dogs$mvo~dogs$lvp)
print(summary(m1.lm))
shapiro.test(rstudent(m1.lm))
qqnorm(rstudent(m1.lm)); qqline(rstudent(m1.lm)); 
hist(rstudent(m1.lm));
m=m1.lm$coefficients  #contient les coef sur toutes les obs
# dogs.fun est la fonction pour calculer les bootsrap
dogs.fun=function(data,i,estim)
{
dt=data[i,]
print(i)
print(dt)
estim=lm(dt$mvo~dt$lvp)$coefficients
c(estim)
}

dogs.boot=boot(dogs,dogs.fun,R=9,stype="i")
print(dogs.boot)
print(dogs.boot$t) #on affiche les estimations ontenues par bootstrap
m1=mean(dogs.boot$t[,1])
m2=mean(dogs.boot$t[,2])
############################

####################### EXERCICE 2

library(MASS)

data(nodal)
nodal$m=factor(nodal$m)
nodal$r=factor(nodal$r)
nodal$aged=factor(nodal$aged)
nodal$stage=factor(nodal$stage)
nodal$grade=factor(nodal$grade)
nodal$xray=factor(nodal$xray)
nodal$acid=factor(nodal$acid)
attach(nodal)
nodal.lm=glm(r~aged+stage+grade+xray+acid,family=binomial)
anova(nodal.lm,test="Chisq")

# nodal.fun est la fonction pour calculer les bootsrap
nodal.fun=function(data,i,estim)
{
dt=data[i,]
print(i)
#print(dt)
estim=glm(dt$r~dt$aged+dt$stage+dt$grade+dt$xray+dt$acid,family=binomial)$coefficients
c(estim)
}

nodal.boot=boot(nodal,nodal.fun,R=99,stype="i")
print(nodal.boot)
print(nodal.boot$t) #on affiche les estimations ontenues par bootstrap


################### EXERCICE 3

######################"
air.fun = function(data)
{    ybar = mean(data$hours)
     para = c(log(ybar),mean(log(data$hours)))
     ll = function(k) {
          if (k <= 0) out <- 1e200 # not NA
          else out <- lgamma(k)-k*(log(k)-1-para[1]+para[2])
         out
     }
     khat = nlm(ll,ybar^2/var(data$hours))$estimate
     c(ybar, khat)
}

air.rg = function(data, mle){    out <- data
     out$hours <- rexp(nrow(out), 1/mle)
     out
}
air.boot = boot(aircondit, air.fun, R=999, sim="parametric",
     ran.gen=air.rg, mle=mean(aircondit$hours))

# The bootstrap p-value can then be approximated by
sum(abs(air.boot$t[,2]-1) > abs(air.boot$t0[2]-1))/(1+air.boot$R)