The French version of this document is no longer maintained: be sure to check the more up-to-date English version.
Utilité des simulations : bootstrap paramétrique
Bootstrap non paramétrique
Exemples particuliers
A FAIRE : effacer ce qui suit ?
A FAIRE
On qualifie de méthode de Monté Carlo tout algorithme de calcul reposant sur des tirages aléatoires, dont le résultat change selon ces tirage, et dont le résultat est d'autant plus précis que le nombre de tirages est élevé.
Par exemple, on peut calculer une valeur approchée de Pi de la manière suivante : on prend des points au hasard (uniformément) dans le carré [0,1]x[0,1] et on regarde combien sont dans le quart de disque de centre O et de rayon 1. Le rapport du nombre de points à l'intérieur sur le nombre de points total est proche de Pi/4.
pi.monte.carlo.plot <- function (n=10000, pch='.', ...) {
x <- runif(n)
y <- runif(n)
interieur <- x^2 + y^2 <= 1
p <- 4*sum(interieur)/n
xc <- seq(0,1,length=200)
yc <- sqrt(1-xc^2)
plot( xc, yc, type='l' )
lines( c(0,1,1,0,0), c(0,0,1,1,0) )
abline(h=0, lty=3)
abline(v=0, lty=3)
points(x[interieur], y[interieur], col='red', pch=pch, ...)
points(x[!interieur], y[!interieur], pch=pch, ...)
title(main=paste("Simulation de Monté Carlo : pi=",p,sep=''))
}
pi.monte.carlo.plot(100, pch='+', cex=3)
pi.monte.carlo.plot()
On peut l'utiliser, comme ici, pour calculer des intégrales multiples. En petites dimensions, ça n'est pas très utile (c'est trop lent), par contre, en dimension supérieures, c'est plus intéressant. Par exemple, pour calculer de manière approchée l'intégrale d'une fonction de 100 variables sur l'hypercube [0,1]^100, la méthode des rectangles, calculer la fonction en chacun des sommets de l'hypercube et faire la moyenne, n'est pas réaliste, car il y a 2^100 sommets ; par contre, une simulation de Monté Carlo donnera un résultat acceptable en un temps raisonnable.
Les simulations permettent aussi d'obtenir l'espérance d'une statistique donnée (je rappelle qu'une statistique, c'est une fonction que l'on évalue sur un échantillon (en toute rigueur, une famille de fonctions (fn : R^n --> R)_n), comme la moyenne, le maximum ou l'écart-type). Considérons par exemple l'expérience suivante : on lance un dé ; si on obtient 1, on tire un nombre au hasard selon une distribution normale de moyenne -1 et d'écart-type 1, sinon, on tire un nombre au hasard selon une distribution normale de moyenne 1 de d'écart-type 2 ; on recommence dix fois. Que dire de la moyenne ou des quartiles de ces dix nombres ?
On peut demander à l'ordinateur de faire des simulations :
n <- 10 x <- ifelse( runif(n)>1/6, rnorm(n,1,2), rnorm(n,-1,1) ) mean(x) quantile(x)
Si on fait plusieurs simulations, on obtient une estimation de l'espérance de la moyenne de des écarts-types) :
n <- 10
N <- 1000
r <- matrix(NA, nr=N, nc=6)
for (i in 1:N) {
x <- ifelse( runif(n)>1/6, rnorm(n,1,2), rnorm(n,-1,1) )
r[i,] <- c( mean(x), quantile(x) )
}
colnames(r) <- c("mean", "0% (min)", "25%", "50% (median)", "75%", "100% (max)")
res <- rbind( apply(r,2,mean), apply(r,2,sd) )
rownames(res) <- c("mean", "sd")Il vient :
mean 0% (min) 25% 50% (median) 75% 100% (max)
mean 0.6365460 -2.256153 -0.6387419 0.5499527 1.8228999 3.807241
sd 0.6446216 1.028375 0.7504127 0.7864454 0.8537194 1.236234On peut même demander un peu plus et regarder la distribution de ces différentes statistiques.
n <- 10
N <- 1000
r <- matrix(NA, nr=N, nc=6)
for (i in 1:N) {
x <- ifelse( runif(n)>1/6, rnorm(n,1,2), rnorm(n,-1,1) )
r[i,] <- c( mean(x), quantile(x) )
}
colnames(r) <- c("mean", "0% (min)", "25%", "50% (median)", "75%", "100% (max)")
rr <- apply(r,2,density)
# Récupérer les valeurs min et max de x et y...
xlim <- range( sapply(rr,function(a){range(a$x)}) )
ylim <- range( sapply(rr,function(a){range(a$y)}) )
plot(NA, xlim=xlim, ylim=ylim, ylab='density')
for (i in 1:6) {
lines(rr[[i]], col=i)
}
legend(par('usr')[2], par('usr')[4], xjust=1, yjust=1,
c('moyenne', "min", "1er quartile", "médiane", "3e quartile", "max"),
lwd=1, lty=1,
col=1:6 )
On voit en particulier sur cette simulation que la moyenne et la médiane ont à peu prés la même valeur, mais la variance de la moyenne est plus petite, i.e., la mpyenne est plus précise : face à ce genre de données, on préfèrerait donc la moyenne à la médiane.
Avoir une idée de la distribution de la statistique qu'on étudie permet de repérer certaines pathologies : ainsi, on suppose souvent que les statistiques qu'on étudie ont des distributions normales ou presque. A défaut de démonstration, des simulations permettent de démontrer que c'est bien le cas. Ainsi, dans l'exemple précédent, on pourrait faire des graphiques quantiles-quantiles.
op <- par(mfrow=c(2,3))
for (i in 1:6) {
qqnorm(r[,i], main=colnames(r)[i])
qqline(r[,i], col='red')
text( par("usr")[1], par("usr")[4], adj=c(-.2,2),
round(shapiro.test(r[,i])$p.value, digits=4) )
}
par(op)
Pour éviter d'avoir à écrire toujours le même code, nous utiliserons désormais la fonction suivante pour nos simulations (dans l'exemple suivant, on regarde la distribution des coefficients d'une régression) :
my.simulation <- function (get.sample, statistic, R) {
res <- statistic(get.sample())
r <- matrix(NA, nr=R, nc=length(res))
r[1,] <- res
for (i in 2:R) {
r[i,] <- statistic(get.sample())
}
list(t=r, t0=apply(r,2,mean))
}
r <- my.simulation(
function () {
n <- 200
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- 1 - x1 + 2 * x2 + rnorm(n)
data.frame(y,x1,x2)
},
function (d) {
lm(y~x1+x2, data=d)$coef
},
R=999
)
matdensityplot <- function (r, ylab='density', ...) {
rr <- apply(r,2,density)
n <- length(rr)
# Récupérer les valeurs min et max de x et y...
xlim <- range( sapply(rr,function(a){range(a$x)}) )
ylim <- range( sapply(rr,function(a){range(a$y)}) )
plot(NA, xlim=xlim, ylim=ylim, ylab=ylab, ...)
for (i in 1:n) {
lines(rr[[i]], col=i)
}
}
matdensityplot(r$t)
On peut aussi utiliser ces simulations pour étudier des statistiques un peu plus complexes comme des différences ou des quotients de moyennes.
q <- runif(1,0,10)
m1 <- runif(1,0,10)
m2 <- q*m1
r <- my.simulation(
function () {
n1 <- 200
n2 <- 100
x1 <- m1*rlnorm(n1)
x2 <- m2*rlnorm(n2)
data.frame( x=c(x1,x2), c=factor(c(rep(1,n1),rep(2,n2))))
},
function (d) {
a <- tapply(d[,1],d[,2],mean)
a[2]/a[1]
},
R=999
)
hist(r$t, probability=T, col='light blue',
main="Distribution du quotient de deux moyennes")
lines(density(r$t), col='red', lwd=3)
abline(v=c(q,r$t0), lty=3, lwd=3, col=c('blue','red'))
legend( par("usr")[2], par("usr")[4], xjust=1, yjust=1,
c("Moyenne théorique", "Moyenne expérimentale"),
lwd=1, lty=3, col=c('blue','red') )
Au lieu de se contenter de calculer la moyenne et l'écart-type de notre statistique, on peut calculer d'autres quantités, comme des intervalles de confiance et des p-valeurs.
# Intervalle de confiance à 5% de l'exemple précédent
hist(r$t, probability=T, col='light blue',
main="Distribution du quotient de deux moyennes")
qt <- quantile(r$t, c(.025,.975))
d <- density(r$t)
o <- d$x>qt[1] & d$x<qt[2]
lines(d$x[o], d$y[o], col='red', lwd=3)
lines(d, col='red')
abline(v=c(q,r$t0), lty=3, lwd=3, col=c('blue','red'))
legend( par("usr")[2], par("usr")[4], xjust=1, yjust=1,
c("Moyenne théorique", "Moyenne expérimentale",
"Intervalle de confiance à 5%"),
lwd=c(1,1,3), lty=c(3,3,1), col=c('blue','red', 'red') )
A FAIRE : mettre le code suivant plus loin, dans la partie sur le
bootstrap non paramétrique
# On a deux échantillons, de tailles différentes, et on cherche le
# quotient de leurs moyennes
n1 <- 200
n2 <- 100
q <- runif(1,0,10)
m1 <- runif(1,0,10)
m2 <- q*m1
x1 <- m1*rlnorm(n1)
x2 <- m2*rlnorm(n2)
d <- data.frame( x=c(x1,x2), c=factor(c(rep(1,n1),rep(2,n2))))
R <- 999
r <- boot(d,
function (d,i) {
a <- tapply(d[,1][i],d[,2][i],mean)
a[2]/a[1]
},
R=R)
hist(r$t, probability=T, col='light blue',
main="Distribution du quotient de deux moyennes")
qt <- quantile(r$t, c(.025,.975))
d <- density(r$t)
o <- d$x>qt[1] & d$x<qt[2]
lines(d$x[o], d$y[o], col='red', lwd=3)
lines(d, col='red')
abline(v=c(q,r$t0), lty=3, lwd=3, col=c('blue','red'))
legend( par("usr")[2], par("usr")[4], xjust=1, yjust=1,
c("Moyenne théorique", "Moyenne expérimentale",
"Intervalle de confiance à 5%"),
lwd=c(1,1,3), lty=c(3,3,1), col=c('blue','red', 'red') )
Les simulations ne permettent pas seulement de regarder la distribution des estimations des paramètres du modèle, mais aussi la distribution des prédictions.
A FAIRE : un exemple
Calcul de l'erreur quadratique moyenne pour une régression classique.
A FAIRE : un exemple avec du bootstrap paramétrique, sans la commande boot.
# Les variables à prédire sont les premières du data.frame
my.simulation.predict <- function (
get.training.sample,
get.test.sample=get.training.sample,
get.model,
get.predictions=predict,
R=999) {
r <- matrix(NA, nr=R, nc=3)
colnames(r) <- c("biais", "variance", "erreur quadratique")
for (i in 1:R) {
d.train <- get.training.sample()
d.test <- get.test.sample()
m <- get.model(d.train)
p <- get.predictions(m, d.test)
if(is.vector(p)){ p <- data.frame(p) }
d.test <- d.test[,1:(dim(p)[2])]
b <- apply(d.test-p, 2, mean)
v <- apply(d.test-p, 2, var)
e <- apply((d.test-p)^2, 2, mean)
r[i,] <- c(b,v,e)
}
list(t=r, t0=apply(r,2,mean))
}
get.sample <- function () {
n <- 200
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- sin(x1) + cos(x2) - x1*x2 + rnorm(n)
data.frame(y,x1,x2)
}
r <- my.simulation.predict(
get.sample,
get.model = function (d) {
lm(y~x1+x2,data=d)
}
)
hist(r$t[,1], probability=T, col='light blue',
main="Biais d'une régression linéaire")
lines(density(r$t), col='red', lwd=3)
abline(v=r$t0[1], lty=3)
A FAIRE : un autre exemple pour une méthode biaisée. (Mais il faut quelque chose de vraiment biaisé...)
En estimant l'erreur quadratique moyenne pour plusieurs méthodes, on peut choisir la meilleure.
get.sample <- function () {
n <- 200
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- sin(x1) + cos(x2) - x1*x2 + rnorm(n)
data.frame(y,x1,x2)
}
r1 <- my.simulation.predict(
get.sample,
get.model = function (d) {
lm(y~x1+x2,data=d)
}
)
r2 <- my.simulation.predict(
get.sample,
get.model = function (d) {
lm(y~poly(x1,3)+poly(x2,3),data=d)
}
)
library(splines)
r3 <- my.simulation.predict(
get.sample,
get.model = function (d) {
lm(y~ns(x1)+ns(x2),data=d)
}
)
matdensityplot(cbind(r1$t[,1], r2$t[,1], r3$t[,1]),
main="Biais de trois régressions")
abline(v=c(r1$t0[1],r2$t0[1],r3$t0[1]), lty=3, col=1:3)
legend( par("usr")[2], par("usr")[4], xjust=1, yjust=1,
c("Régression linéaire", "Régression polynomiale", "Splines"),
lty=1, lwd=1, col=1:3)
%--On ne voit pas grand-chose sur le dessin : regardons les erreurs quadratiques moyennes.
> rbind( r1$t0, r2$t0, r3$t0 )
biais variance erreur quadratique
[1,] -0.004572499 2.301057 2.311898
[2,] -0.001475014 2.186536 2.194847
[3,] -0.001067888 2.331862 2.342088La deuxième méthode semble très légèrement préférable : la différence des moyennes est faible mais statistiquement significative.
> t.test( r2$t[,3], r1$t[,3] )
Welch Two Sample t-test
data: r2$t[, 3] and r1$t[, 3]
t = -7.1384, df = 1994.002, p-value = 1.317e-12
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.14920930 -0.08489377
sample estimates:
mean of x mean of y
2.194847 2.311898
> wilcox.test( r2$t[,3], r1$t[,3] )
Wilcoxon rank sum test with continuity correction
data: r2$t[, 3] and r1$t[, 3]
W = 397602, p-value = 3.718e-15
alternative hypothesis: true mu is not equal to 0Voici un autre exemple, non plus pour choisir entre deux ou trois méthodes différentes, mais pour choisir un paramètre dans une méthode.
library(MASS)
get.sample <- function () {
n <- 20
x <- rnorm(n)
x1 <- x + .2*rnorm(n)
x2 <- x + .2*rnorm(n)
x3 <- x + .2*rnorm(n)
x4 <- x + .2*rnorm(n)
y <- x + (x1+x2+x3+x4)/4 + rnorm(n)
data.frame(y,x1,x2,x3,x4)
}
get.model <- function (d) {
r <- lm.ridge(y~., data=d)
}
get.predictions <- function (r,d) {
m <- t(as.matrix(r$coef/r$scales))
inter <- r$ym - m %*% r$xm
drop(inter) + as.matrix(d[,-1]) %*% drop(m)
}
lambda <- c(0,.001,.002,.005,.01,.02,.05,.1,.2,.5,1,2,5,10,20,50,100)
n <- length(lambda)
res <- rep(NA,n)
for (i in 1:n) {
res[i] <- my.simulation.predict(
get.sample,
get.model = function (d) {
lm.ridge(y~., data=d, lambda=lambda[i])
},
get.predictions = get.predictions,
R=99
)$t0[3]
}
plot(res~lambda, type='l', log='x')
abline(h=res[1], lty=3)
i <- (1:n)[ res == min(res) ]
x <- lambda[i]
y <- res[i]
points(x,y,col='red',pch=15)
title(main="Détermination du paramètre d'une ridge regression par boostrap")
Les estimateurs obtenus à l'aide de la méthode du maximum de vraissemblance sont souvent biaisés : l'espérance de l'estimateur n'a pas la bonne valeur.
A FAIRE : un exemple d'estimateur MLE biaisé
Dans d'autres situations ce sont les prédictions qui sont suceptibles d'être biaisées.
On peut aussi se retrouver avec des estimateurs intentionnellement biaisés : dans le cas de la ridge regression ou de la régression sur les composantes principales, on accepte un biais en contrepartie d'une variance plus faible, ce qui fait décroître l'erreur quadratique moyenne.
Des simulations permettent d'estimer le biais d'un estimateur ou d'une prédiction.
Si on reprend l'exemple de la ridge regression, les prédictions ne sont pas vraiment biaisées.
library(MASS)
get.sample <- function () {
n <- 20
x <- rnorm(n)
x1 <- x + .2*rnorm(n)
x2 <- x + .2*rnorm(n)
x3 <- x + .2*rnorm(n)
x4 <- x + .2*rnorm(n)
y <- 1+x1+2*x2+x3+2*x4 + rnorm(n)
data.frame(y,x1,x2,x3,x4)
}
get.model <- function (d) {
r <- lm.ridge(y~., data=d)
}
get.predictions <- function (r,d) {
m <- t(as.matrix(r$coef/r$scales))
inter <- r$ym - m %*% r$xm
drop(inter) + as.matrix(d[,-1]) %*% drop(m)
}
lambda <- c(0,.001,.002,.005,.01,.02,.05,.1,.2,.5,1,2,5,10,20,50,100)
lambda <- exp(seq(-7,7,length=100))
n <- length(lambda)
res <- rep(NA,n)
for (i in 1:n) {
res[i] <- my.simulation.predict(
get.sample,
get.model = function (d) {
lm.ridge(y~., data=d, lambda=lambda[i])
},
get.predictions = get.predictions,
R=20
)$t0[1]
}
plot(res~lambda, type='l', log='x')
abline(h=0,lty=3)
title("Biais des prédictions de la ridge regression")
Par contre, les estimations des coefficients sont biaisées.
get.sample <- function () {
n <- 20
x <- rnorm(n)
x1 <- x + .2*rnorm(n)
x2 <- x + .2*rnorm(n)
x3 <- x + .2*rnorm(n)
x4 <- x + .2*rnorm(n)
y <- 1+2*x1+3*x2+4*x3+5*x4 + rnorm(n)
data.frame(y,x1,x2,x3,x4)
}
s <- function(d,l) {
r <- lm.ridge(y~., data=d, lambda=l)
m <- t(as.matrix(r$coef/r$scales))
inter <- r$ym - m %*% r$xm
c(inter,m)
}
lambda <- exp(seq(-7,7,length=100))
n <- length(lambda)
res <- matrix(NA,nr=n,nc=5)
for (i in 1:n) {
res[i,] <- my.simulation(
get.sample,
function (d) { s(d,lambda[i]) },
R=20
)$t0
}
matplot(lambda,res,type='l',lty=1,log='x')
abline(h=1:5,lty=3,col=1:5)
title("Biais des coefficients de la ridge regression")
Prenons un autre exemple, la régression sur les composantes principales, et regardons la valeur du biais en fonction de la valeur des paramètres.
get.sample <- function (a,b) {
n <- 20
x <- rnorm(n)
x1 <- x + .2* rnorm(n)
x2 <- x + .2* rnorm(n)
y <- a*x1 + b*x2 + rnorm(n)
data.frame(y,x1,x2)
}
get.parameters <- function (d) {
y <- d[,1]
x <- d[,-1]
p <- princomp(x)
r <- lm(y ~ p$scores[,1] -1)
drop(p$loadings %*% c(r$coef,0))
}
N <- 5
a.min <- 0
a.max <- 10
res <- matrix(NA,nr=N,nc=N)
a <- seq(a.min,a.max,length=N)
b <- seq(a.min,a.max,length=N)
for (i in 1:N) {
for(j in 1:N) {
res[i,j] <- my.simulation(
function () { get.sample(a[i],b[j]) },
get.parameters,
R=20
)$t0[1] - a[i]
}
}
persp(res)
image(a,b,res,col=topo.colors(255))
text( matrix(a,nr=N,nc=N),
matrix(b,nr=N,nc=N, byrow=T),
round(res) )
On peut utiliser cette connaissance du biais pour tenter de le corriger. Ainsi, dans l'exemple précédent, le biais sur a est à peu près (b-a)/2.
get.parameters <- function (d) {
y <- d[,1]
x <- d[,-1]
p <- princomp(x)
r <- lm(y ~ p$scores[,1] -1)
res <- drop(p$loadings %*% c(r$coef,0))
a <- res[1]
b <- res[2]
c( a - (b-a)/2, b - (a-b)/2 )
}
N <- 5
a.min <- 0
a.max <- 10
res <- matrix(NA,nr=N,nc=N)
a <- seq(a.min,a.max,length=N)
b <- seq(a.min,a.max,length=N)
for (i in 1:N) {
for(j in 1:N) {
res[i,j] <- my.simulation(
function () { get.sample(a[i],b[j]) },
get.parameters,
R=20
)$t0[1] - a[i]
}
}
persp(res)
image(a,b,res,col=topo.colors(255))
text( matrix(a,nr=N,nc=N),
matrix(b,nr=N,nc=N, byrow=T),
round(res) )
Ca ne marche pas du tout : on a une estimation du biais en fonction des valeurs théoriques des paramètres, alors qu'il faudrait une estimation du biais en fonction des valeurs expérimentales.
get.parameters <- function (d) {
y <- d[,1]
x <- d[,-1]
p <- princomp(x)
r <- lm(y ~ p$scores[,1] -1)
drop(p$loadings %*% c(r$coef,0))
}
N <- 5
a.min <- 0
a.max <- 10
res <- matrix(NA,nr=N,nc=N)
res.a <- matrix(NA,nr=N,nc=N)
res.b <- matrix(NA,nr=N,nc=N)
a <- seq(a.min,a.max,length=N)
b <- seq(a.min,a.max,length=N)
for (i in 1:N) {
for(j in 1:N) {
r <- my.simulation(
function () { get.sample(a[i],b[j]) },
get.parameters,
R=20
)$t0
res.a[i,j] <- r[1]
res.b[i,j] <- r[2]
res[i,j] <- r[1] - a[i]
}
}
plot(res.a, res.b, type='n')
text(res.a, res.b, round(res))
plot(as.vector(res) ~ as.vector(res.a-res.b))
Eh bien, ça n'était pas un bon exemple (mais dans certains cas, comme la ridge regression plus haut, on peut avoir une estimation un peu plus fidèle du biais en fonction des estimations des paramètres qui permet de les corriger).
Dans les discussions précédentes, nous avons supposé qu'on disposait d'un modèle décrivant fidèlement les observations. Dans ce cas, on peut généralement calculer exactement l'espérance, la variance, la distribution des estimateurs, leur biais et l'erreur quadratique moyenne des prédictions.
Néanmoins, ces simulations, que l'on qualifie de "bootstrap paramétrique", sont quand même utiles, en particulier quand le modèle est trop complexe pour qu'on puisse faire des calculs exacts ou quand le modèle nous est donné comme une boite noire (on a un programme pour générer de nouvelles observations, mais on ne sait pas comment il fonctionne).
Nous allons voir que des simulations permettent de répondre aux même questions en l'absence de modèle : c'est le bootstrap non paramétrique.
Pour faire des simulations sans modèle, avec juste un échantillon, on n'a pas d'autre choix que d'utiliser cet échantillon.
Une première idée consiste à diviser l'échantillon en deux, en utilisant 75% des données pour trouver l'estimateur ou les paramètres nécessaires aux prédictions, et en utilisant les 25% restants pour vérifier la pertinence des prédictions. On peut ensuite recommencer avec une autre partition des données en deux, de manière à avoir une idée plus fine de la répartition des erreurs de prédictions. Ensuite, si le modèle ou l'algorithme semble convenir, on le relance sur toutes les données (pour avoir un résultat plus précis). C'est ce qu'on appelle la validation croisée.
my.simulation.cross.validation <- function (
s, # sample
k, # number of observations to fit the the model
get.model,
get.predictions=predict,
R = 999 )
{
n <- dim(s)[1]
r <- matrix(NA, nr=R, nc=3)
colnames(r) <- c("biais", "variance", "erreur quadratique")
for (i in 1:R) {
j <- sample(1:n, k)
d.train <- s[j,]
d.test <- s[-j,]
m <- get.model(d.train)
p <- get.predictions(m, d.test)
if(is.vector(p)){ p <- data.frame(p) }
d.test <- d.test[,1:(dim(p)[2])]
b <- apply(d.test-p, 2, mean)
v <- apply(d.test-p, 2, var)
e <- apply((d.test-p)^2, 2, mean)
r[i,] <- c(b,v,e)
}
list(t=r, t0=apply(r,2,mean))
}
n <- 200
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
y <- 2 - x1 - x2 - x3 + rnorm(n)
s <- data.frame(y,x1,x2,x3)
r <- my.simulation.cross.validation(s, 150, function(s){lm(y~.,data=s)}, R=99)
hist( r$t[,3], probability=T, col='light blue',
main="Estimation de l'erreur quadratique moyenne par validation croisée" )
lines(density(r$t[,3]),col='red',lwd=3)
abline(v=mean(r$t[,3]),lty=2,lwd=3,col='red')
La partition des observations n'est pas nécessairement 75%/25%, on peut prendre 90%/10% ou même n'enlever qu'une seule observation : ce dernier cas constitue le Jackknife.
my.simulation.jackknife <- function (
s,
get.model,
get.predictions=predict,
R=999)
{
n <- dim(s)[1]
if(R<n) {
j <- sample(1:n, R)
} else {
R <- n
j <- 1:n
}
p <- get.predictions(get.model(s), s)
if( is.vector(p) ) {
p <- as.data.frame(p)
}
r <- matrix(NA, nr=R, nc=dim(s)[2]+dim(p)[2])
colnames(r) <- c(colnames(s), colnames(p))
for (i in j) {
d.train <- s[-i,]
d.test <- s[i,]
m <- get.model(d.train)
p <- get.predictions(m, d.test)
r[i,] <- as.matrix(cbind(s[i,], p))
}
r
}
n <- 200
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
y <- 2 - x1 - x2 - x3 + rnorm(n)
s <- data.frame(y,x1,x2,x3)
r <- my.simulation.jackknife(s, function(s){lm(y~.,data=s)})
hist( r[,5]-r[,1], col='light blue' )
A FAIRE : Expiquer comment abaisser le biais d'un estimateur à l'aide du Jackknife Expliquer le principe, donner la formule qui en découle, en déduire la formule effectivement utilisée pour le calcul. Donner un exemple (probablement la variance de la population, qui est un estimateur biaisé : le bootstrap redonnera la variance de l'échantillon).
Dans la validation croisée (ou le jackknive), chaque observation était utilisée exactement une fois : on construit l'échantillon sur lequel on va faire les calculs à l'aide d'un tirage SANS remise. On peut utiliser un tirage AVEC remise : c'est le principe du bootstrap. Cela revient à faire jouer à l'échantillon le rôle de la population.
On prend un tirage AVEC remise afin que les variables aléatoires obtenues soient indépendantes. (Dans un tirage sans remise, après avoir tiré au hasard une première observation, la seconde qu'on tire au hasard dépend de la première, car on sait qu'elle est différente de la première.)
Pour programmer cela, on pourrait reprendre la fonction que nous avons utilisée pour la validation croisée.
my.simulation.bootstrap <- function (
s,
get.model,
get.predictions=predict,
R = 999 )
{
n <- dim(s)[1]
r <- matrix(NA, nr=R, nc=3)
colnames(r) <- c("biais", "variance", "erreur quadratique")
for (i in 1:R) {
j <- sample(1:n, n, replace=T)
d.train <- s[j,]
d.test <- s
m <- get.model(d.train)
p <- get.predictions(m, d.test)
if(is.vector(p)){ p <- data.frame(p) }
d.test <- d.test[,1:(dim(p)[2])]
b <- apply(d.test-p, 2, mean)
v <- apply(d.test-p, 2, var)
e <- apply((d.test-p)^2, 2, mean)
r[i,] <- c(b,v,e)
}
list(t=r, t0=apply(r,2,mean))
}
n <- 200
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
y <- 2 - x1 - x2 - x3 + rnorm(n)
s <- data.frame(y,x1,x2,x3)
r <- my.simulation.bootstrap(s, function(s){lm(y~.,data=s)}, R=99)
hist( r$t[,3], col='light blue',
main="Estimation de l'erreur quadratique moyenne par bootstrap" )
lines(density(r$t[,3]),col='red',lwd=3)
abline(v=mean(r$t[,3]),lty=2,lwd=3,col='red')
Mais en fait, il y a déjà une fonction qui fait ça.
Voici par exemple notre première utilisation du bootstrap : regarder la distribution d'une statistique (ici, la moyenne d'un échantillon).
library(boot)
n <- 200
x <- rlnorm(n)
r <- boot(x, function(s,i){mean(s[i])}, R=99)
hist(r$t, probability=T, col='light blue')
lines(density(r$t),col='red',lwd=3)
abline(v=mean(r$t),lty=2,lwd=3,col='red')
title("Estimation par bootstrap de la distribution de la moyenne d'un échantillon")
On donne les observations en premier argument de la fonction boot, sous forme d'un vecteur (s'il n'y a qu'une seule variable), d'un tableau ou d'un data.frame. Le second argument est la statistique, fonction qui prend les données et un vecteur d'indices et renvoie un vecteur. Le dernier argument est le nombre d'échantillons de bootstrap.
n <- 200
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
y <- -2 - x1 +0* x2 + x3 + rnorm(n)
s <- data.frame(y,x1,x2,x3)
r <- boot(s, function(s,i){lm(y~.,data=s[i,])$coef}, R=99)
plot(NA,xlim=range(r$t), ylim=c(0,6.5))
for (i in 1:4) {
lines(density(r$t[,i]), lwd=2, col=i)
}
Le résultat est une liste dont essentiellement deux composantes nous intéresseront : t, qui contient toutes les estimations de la statisque à étudier et t0 qui contient l'estimation pour tout l'échantillon. La fonction print.boot donne aussi le biais :
> r
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = s, statistic = function(s, i) {
lm(y ~ ., data = s[i, ])$coef
}, R = 99)
Bootstrap Statistics :
original bias std. error
t1* -2.0565761 -0.0001888077 0.06901359
t2* -1.0524309 -0.0060001140 0.06248444
t3* -0.1340356 -0.0082527142 0.06299820
t4* 1.0176907 0.0013180553 0.06206947Regardons maintenant l'autre utilisation du bootstrap : l'analyse des prédictions. Il s'agit simplement d'étudier la statistique "erreur quadratique moyenne des prévisions".
n <- 200
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
y <- 1+x1+x2+x3+rnorm(n)
d <- data.frame(y,x1,x2,x3)
r <- boot(d,
function (d,i) {
r <- lm(y~.,data=d[i,])
p <- predict(r, d[-i,])
mean((p-d[,1][-i])^2)
},
R=99)
hist(r$t,probability=T,col='light blue',
main="Répartition de l'erreur quadratique moyenne")
lines(density(r$t),lwd=3,col='red')
abline(v=mean(r$t),lty=3)
Le crédo sous-jacent à cette méthode est le suivant : les liens entre l'échantillon dont on dispose et la population sont à peu près les mêmes que ceux entre cet échantillon et des échantillons de même taille construits par tirage avec remise.
Vérifons sur une simulation que cette hypothèse n'est pas trop farfelue.
get.sample <- function (n=50) { rnorm(n) }
get.statistic <- function (x) { mean(x) }
# Traçons la densité de l'estimateur
R <- 200
x <- rep(NA,R)
for (i in 1:R) {
x[i] <- get.statistic(get.sample())
}
plot(density(x), ylim=c(0,4), lty=1,lwd=3,col='red',
main="Pertinence du bootstrap")
# Traçons quelques estimations de cette densité obtenus par le bootstrap
for (i in 1:5) {
x <- get.sample()
r <- boot(x, function(s,i){mean(s[i])}, R=R)
lines(density(r$t))
}
curve(dnorm(x,sd=1/sqrt(50)),add=T, lty=2,lwd=3,col='blue')
legend(par('usr')[2], par('usr')[4], xjust=1, yjust=1,
c('bootstrap paramétrique', 'bootstraps non paramétriques',
"Courbe théorique"),
lwd=c(3,1,2),
lty=c(1,1,2),
col=c("red", par("fg"), "blue"))
Moralité : le bootstrap ne fait pas de miracle.
Le bootstrap n'est pertinent que pour certains estimateurs, pas trop pathologiques : en particulier, il ne faut pas prendre des choses comme le maximum, sous peine d'avoir des résultats complètement farfelus.
n <- 50
x <- rnorm(n)
r <- boot(x, function(s,i){max(s[i])}, R=999)
hist(r$t, probability=T, col='light blue',
main="Une statistique pour laquelle le bootstrap n'est pas très pertinent...")
lines(density(r$t,bw=.05), col='red', lwd=3)
La commande boot.ci permet d'obtenir des intervalles de confiance.
n <- 20
x <- rnorm(n)
r <- boot(x, function(x,i){mean(x[i])}, R=99)
> boot.ci(r)
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 99 bootstrap replicates
CALL :
boot.ci(boot.out = r)
Intervals :
Level Normal Basic
95% (-0.1457, 0.6750 ) (-0.2429, 0.5969 )
Level Percentile BCa
95% (-0.1050, 0.7348 ) (-0.0559, 0.7891 )
Calculations and Intervals on Original Scale
Some basic intervals may be unstable
Some percentile intervals may be unstable
Warning : BCa Intervals used Extreme Quantiles
Some BCa intervals may be unstable
Warning messages:
1: Bootstrap variances needed for studentized intervals in: boot.ci(r)
2: Extreme Order Statistics used as Endpoints in: norm.inter(t, adj.alpha)
n <- 20
N <- 100
r <- matrix(NA,nr=N,nc=2)
for (i in 1:N) {
x <- rnorm(n)
r[i,] <- boot.ci(boot(x,function(x,i){mean(x[i])},R=99))$basic[c(4,5)]
}
plot(NA,xlim=c(-1.5,1.5), ylim=c(0,2),
main="Répartition des bornes des intervalles de confiance donnés par le bootstrap")
lines(density(r[,1]), lwd=3, col='red')
lines(density(r[,2]), lwd=3, col='blue')
abline(v=0,lty=3)
Le bootstrap est biaisé (son estimation du biais est biaisée vers zéro), car certaines observations sont utilisées à la fois dans l'échantillon pour construire le modèle et dans l'échantillon pour le valider. Le bootstrap "hors du sac" (out-of-the-bag) et le bootstrap .632 tentent de corriger ce biais.
Le bootstrap "hors du sac" consiste à ne pas utiliser toutes les observations pour valider le modèle mais uniquement celles qui ne figurent pas déjà dans l'échantillon ayant servi à le construire (c'est d'ailleurs ce qu'on faisait pour la validation croisée).
En fait, le bootstrap "out-of-the-bag" est quand-même biaisé, mais dans l'autre sens. Pour tenter de corriger ce biais, on peut faire une moyenne pondérée du bootstrap initial et du bootstrap oob.
.368 * (biais estimé par le bootstrap) + .632 * (biais estimé par le bootstrap oob)
(le coefficient .632 s'interprète ainsi : pour n grand, les échantillons de bootstrap contiennent en moyenne 63,2% des observations initiales).
A FAIRE : exemple ?
L'exemple qui suit n'est pas bon : l'estimateur n'est pas biaisé...
n <- 200
x1 <- rnorm(n)
x2 <- rnorm(n)
x3 <- rnorm(n)
y <- 1+x1+x2+x3+rnorm(n)
d <- data.frame(y,x1,x2,x3)
r1 <- boot(d, function (d,i) {
r <- lm(y~.,data=d[i,])
p <- predict(r, d)
mean(p-d[,1])
},
R=99)$t
r2 <- boot(d, function (d,i) {
r <- lm(y~.,data=d[i,])
p <- predict(r, d[-i,])
mean(p-d[,1][-i])
},
R=99)$t
r3 <- .368*r1 + .632*r2
hist(r1,probability=T,col='light blue',
main="Estimations du biais : bootstrap, oobb, .632")
lines(density(r1),lwd=3,col='red')
lines(density(r2),lwd=3,col='green')
lines(density(r3),lwd=3,col='blue')
legend(par('usr')[2], par('usr')[4], xjust=1, yjust=1,
c("bootstrap classique", "bootstrap hors-du-sac", "boootstrap .632"),
lwd=3,
col=c("red", "green", "blue")
)
Face à une régression linéaire, on a tendance à bootstrapper les résidus.
A FAIRE : comparaison des trois méthodes bootstrap paramétrique bootstrap normal bootstrap des résidus pour un modèle Y = a X1 + b X2
Le bootstrap ne s'applique pas directement aux séries temporelles, car leurs observations ne sont pas indépendantes. Généralement, on tire au hasard des morceaux de série temporelle (tirage avec remise) et on les met bout-à-bout.
A FAIRE : dans les exemples ci-dessus, il faut utiliser la fonction boot, même pour le bootstrap paramétrique...
get.sample.ridge <- function (n=100) {
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- 1 - x1 - x2 + rnorm(n)
data.frame(y,x1,x2)
}
statistic.ridge <- function (x, lambda=0) {
r <- lm.ridge(y~x1+x2, data=x, lambda=lambda)
}
statistic.ridge.for.boot <- function (x,i, lambda=0) {
statistic.ridge(x[i,], lambda=10)
}
# Bootstrap non paramétrique
n <- 20
N <- 100
x <- get.sample.ridge()
r <- boot(x, statistic.ridge.for.boot, R=999)Il vient :
> r
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = x, statistic = statistic.ridge.for.boot, R = 999)
Bootstrap Statistics :
original bias std. error
t1* 1.1380068 0.0007701467 0.10555253
t2* -1.0200495 0.0037754155 0.09702321
t3* -0.9293355 0.0017206597 0.08876848Et pour le bootstrap paramétrique :
# Bootstrap paramétrique
ran.ridge <- function (d, p=c(1,-1,-1)) {
n <- 100
x1 <- rnorm(n)
x2 <- rnorm(n)
y <- p[0] + p[1]*x1 - p[2]*x2 + rnorm(n)
data.frame(y,x1,x2)
}
r <- boot(x, statistic.ridge.for.boot, R=999, sim='parametric',
ran.gen=ran.ridge, mle=c(1,-1,-1))
A FAIRE : il y a une erreur...
A FAIRE : estimation du biais pour des prédictions biaisées (prendre une régression normale, faire les prédictions et ajouter 1 -- comme ça, se sera biaisé). Utiliser la commande boot pour faire ces simulations... Tracer la courbe représentant la densité de l'erreur de prédiction pour la méthode choisie (régression intentionnelement biaisée), et pour cette même méthode après correction du biais.
A FAIRE : comparaison de plusieurs régressions.
y <- sin(2*pi*x)+.3*rnorm(n) y~x y~ns(x) y~poly(x,20) Comparaison : estimation de l'erreur quadratique par bootstrap, AIC et BIC
Diviser les observations en dix groupes : 1, 2, 3,..., 10. Utiliser les groupes 2,3,...,10 pour calculer l'estimateur, vérifier avec le groupe 1, en déduire une estimation du biais. Recommencer avec les groupes 1,3,4,...,10 (i.e., en mettant le groupe 2 de côté), etc. En faisant la moyenne de ces dix biais, on obtient une estimation du biais. Calculer l'estimateur avec toutes les données, le corriger à l'aide du biais précédent.
A FAIRE : réutiliser le paragraphe suivant.
Nous nous plaçons dans la situation suivante : on a un échantillon, d'une population non normale, et on cherche à estimer l'un des paramètres de cette population (par exemple, la moyenne, ou le maximum, ou les coefficients d'une régression). Quel estimateur choisir ? Est-il biaisé ? Peut-on avoir une estimation de ce biais ? Quelle est sa variance ? Suit-il une distribution normale ? Peut-on avoir des intervalles de confiance ? Sans rien savoir sur la population, on ne peut pas dire grand-chose. Mais on sait quelque-chose sur la population : on en a un échantillon. L'idée est d'utiliser cet échantillon pour construire d'autres échantillons (on fait simplement des tirages AVEC remise -- c'est la grosse différence avec la validation croisée) : on peut alors regarder le comportement de notre estimateur sur ces échantillons. Ca donne une idée du comportement de l'estimateur.
A FAIRE : dire qu'on ne peut pas étudier n'importe quelle statistique à l'aide du bootstrap : il faut qu'elle possèque certaines propriétés théoriques.
A FAIRE : Réécrire complètement ce qui suit. Il faut prendre un exemple autre que le maximum : toutes les fonctions statistiques ne se prètent pas au bootstrap, il faut qu'elles vérifient certaines propriétés théoriques. Or le maximum est l'exemple typique de la fonction statistique qui ne marche pas.
A FAIRE : biais et optimisme
Expliquons la différence entre le biais et l'optimisme. On dispose d'une population (que l'on ne connait pas), d'un échantillon, et de plusieurs échantillons de bootstrap. On dispose aussi d'un estimateur, sensé mesurer un certain paramètre de la population. Le biais de l'estimateur est l'espérance de la différence entre la valeur de ce paramètre sur la population et la valeur de l'estimateur sur l'échantillon ; l'optimisme est la différence entre la valeur de l'estimateur sur l'échantillon et la moyenne des valeurs de cet estimateur sur les échantillons de bootstrap. Le biais est une comparaison entre l'estimateur et le paramètre qu'il est sensé mesurer, alors que l'optimisme ne fait intervenir que l'estimateur, pas le paramètre.
A FAIRE : exercice : écrire une fonction qui fait du bootstrap.
my.bootstrap <- function (data, estimator, size, ...) {
values <- NULL
n <- dim(data)[1]
for (i in 1:size) {
bootstrap.sample <- data[ sample(1:n, n, replace=T), ]
values <- rbind(values, t(estimator(bootstrap.sample, ...)))
}
values
}A FAIRE : mentionner "bootstrap" :
(Il y a aussi une bibliothèque "bootstrap", mais elle est beaucoup plus vieille.)
A FAIRE : présentation de la commande boot.
La commande boot attend en argument : les données (généralement, un vecteur ou un data.frame), l'estimateur (une fonction qui prend en arguments les données et les indices des observations de l'échantillon de bootstrap) et le nombre d'échantillons de bootstrap à calculer.
library(boot)
my.max <- function (x, i) { max(x[i]) }
r <- boot(faithful$eruptions, my.max, R=100)
plot(r)
> r
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = faithful$eruptions, statistic = my.max, R = 100)
Bootstrap Statistics :
original bias std. error
t1* 5.1 -0.02265 0.0367342
Je ne sais pas ce que c'est... (C'est pour regarder l'influence de chacune des observations sur le résultat.)
A FAIRE
jack.after.boot(r)
L'efficacité du bootstrap repose sur la croyance que les liens entre la population (qu'on ne connait pas) et les échantillons (on en a un) sont les mêmes que ceux entre l'échantillon dont on dispose et ses échantillons de bootstrap. A priori, le lien est très faible.
Pour tenter de justifier le bootstrap, on peut faire l'expérience suivante.
1. Choisir une distribution 2. Prendre un échantillon 3. Appliquer le bootstrap pour estimer un paramètre ou des intervalles de confiance 4. Recommencer en 2.
Voici ce que cela peut donner :
A FAIRE : tentative de justification du bootstrap
On peut aussi utiliser ça quand on dispose de données réelles, pour lesquelles on a un modèle : on cherche les paramètres de ce modèle, par exemple, à l'aide de la méthode du maximum de vraissemblance, et on prend des échantillons à l'aide de ce modèle et non plus par tirages avec remise.
A FAIRE : exemple de bootstrap paramétrique (reprendre la durée des éruptions du geyser, qu'on avait modélisée comme un mélange de deux gaussiennes)
A FAIRE : comparaison du bootstrap paramétrique et non paramétrique. (en particulier, regarder combien d'échantillons il faut pour avoir des estimateurs fiables).
Il y a un gros problème dans le principe même du bootstrap : il y a des observations qui sont utilisées à la fois lors de la modélisation et lors de la validation. Cela risque de donner une estimation du biais biaisée vers zéro. Pour y remédier, on peut valider une modélisation effectuée sur un échantillon de bootstrap uniquement sur les observations qui ne sont pas dedans : c'est le "out-of-the-bag bootstrap" (A FAIRE : traduire).
A FAIRE : relire l'exemple suivant
get.sample <- function (n=50) {
x <- rnorm(n)
y <- 1-x+rnorm(n)
data.frame(x,y)
}
s <- get.sample()
R <- 100
n <- dim(s)[1]
x <- s$x
y <- s$y
res <- rep(NA,R)
res.oob <- rep(NA,R)
for (k in 1:R) {
i <- sample(n,n,replace=T)
xx <- x[i]
yy <- y[i]
r <- lm(yy~xx)
ip <- rep(T,n)
ip[i] <- F
xp <- x[ip]
yp <- predict(r,data.frame(xx=xp))
try( res.oob[k] <- sd(yp - s$y[ip]) )
res[k] <- sd( yy - predict(r) )
}
plot(density(res), main="Basic and out-of-the-box bootstrap")
lines(density(res.oob),lty=2)
abline(v=1,lty=3)
legend( par("usr")[2], par("usr")[4], xjust=1, yjust=1,
c("Basic bootstrap", "Out-of-the-box bootstrap"),
lwd=1, lty=c(1,2) )
On obtient :
> mean(res) [1] 0.9221029 > mean(res.oob) [1] 0.9617512
au lieu de 1. Dans bien des cas, cependant, le bootstrap classique donne des estimations trop optimistes alors que l'autre est trop pessimiste : on fait souvent une moyenne des deux (le coefficient .632 s'interprète ainsi : pour n grand, les échantillons de bootstrap contiennent en moyenne 63,2% des observations initiales).
> .368*mean(res) + .632*mean(res.oob) [1] 0.9471606
A FAIRE : relire ce qui suit
A FAIRE : changer le titre de cette section ?
Nous venons de voir comment le bootstrap permettait d'estimer le biais d'un estimateur. Il permet aussi de comparer différents estimateurs afin de choisir celui qui a la variance la plus faible.
C'est par exemple la situation que nous avons rencontrée quand nous avons parlé de ridge regression : en cas de multicolinéarité, on ne fait plus la régression en utilisant
b = (X'X)^-1 X' y
(car la matrice X'X est presque singulière) mais en utilisant
b = (X'X + kI)^-1 X' y.
Essayons de trouver la valeur de k minimisant la variance de l'estimateur. (Dans cette situation, comme il n'y a rien d'autre que de l'algèbre linéaire et des distributions normales, on doit pouvoir calculer la variance sans passer par le boostrap.)
echantillon.multicolineaire <- function (n=100) {
x0 <- rep(1,n)
x1 <- x0 + .1*rnorm(n)
x2 <- x0 + .1*rnorm(n)
x3 <- x0 + .1*rnorm(n)
y <- x0 + x1 + x2 + x3 + rnorm(n)
x <- cbind(x0,x1,x2,x3)
m <- cbind(y,x)
m
}
m <- echantillon.multicolineaire()
f <- function (m, k) {
x <- m[,-1]
y <- m[,1]
n <- dim(x)[2]
solve( t(x) %*% x + k*diag(rep(1,n)), t(x) %*% y )
}
N <- 200
variance <- rep(NA, N)
biais <- rep(NA, N)
exact <- f(m,0)
k <- exp(seq(.01,4, length=N))-1
for (j in 1:N) {
g <- function (m, i) { f(m[i,], k[j]) }
s <- boot(m, g, R=1000)
variance[j] <- max(var(s$t))
biais[j] <- max(exact - apply(s$t,2,mean))
}
plot(variance ~ k, log='x', type='l')
par(new=T)
plot(biais ~ k, type='l', log='x', col='red', xlab='', ylab='', axes=F)
axis(4, col="red")
legend( par("usr")[2], par("usr")[4], yjust=1, xjust=0,
c("variance", "optimisme"),
lwd=1, lty=1,
col=c(par('fg'), 'red') )
A FAIRE : ce graphique ne suggère pas k=10, comme j'aimerais... Avant d'utiliser la commande boot(), j'obtenais ça :
my.bootstrap <- function (data, estimator, size, ...) {
values <- NULL
n <- dim(data)[1]
for (i in 1:size) {
bootstrap.sample <- data[ sample(1:n, n, replace=T), ]
values <- rbind(values, t(estimator(bootstrap.sample, ...)))
}
values
}
n <- 100
x0 <- rep(1,n)
x1 <- x0 + .1*rnorm(n)
x2 <- x0 + .1*rnorm(n)
x3 <- x0 + .1*rnorm(n)
y <- x0 + x1 + x2 + x3 + rnorm(n)
x <- cbind(x0,x1,x2,x3)
f <- function (m, k) {
x <- m[,-1]
y <- m[,1]
n <- dim(x)[2]
solve( t(x) %*% x + k*diag(rep(1,n)), t(x) %*% y )
}
N <- 200
variance <- rep(NA, N)
biais <- rep(NA, N)
k <- exp(seq(0,6, length=N))-1
for (i in 1:N) {
s <- my.bootstrap(cbind(y,x), f, 100, k[i])
variance[i] <- max(apply(s, 2, var))
biais[i] <- max(abs( s -
matrix( as.vector(f(cbind(y,x),0)), nr=dim(s)[1], nc=dim(s)[2], byrow=T )
))
}
plot(variance ~ k, log='x', type='l')
par(new=T)
plot(biais ~ k, type='l', log='x', col='red', xlab='', ylab='', axes=F)
axis(4)
On choisit k=0.2
> g <- function (m, i) { f(m[i,], .2) }
> s <- boot(m, g, R=100)
> apply(s$t,2,mean)
[1] 1.1294813 1.0082432 0.3821003 1.2957629
> g <- function (m, i) { f(m[i,], 10) }
> s <- boot(m, g, R=100)
> apply(s$t,2,mean)
[1] 0.9402689 0.9762665 0.8696934 0.9346982Ces valeurs sont bien mieux que celles obtenues pour k=0 (i.e., par la régression usuelle) : les coefficients de x1, x2 et x3 sont sensés être égaux.
> as.vector(exact) [1] 1.3112394 1.0350861 0.2463275 1.2259999
On peut vérifier la précision des prédictions.
s <- my.bootstrap(cbind(y,x), f, 100, 10)
b <- apply(s, 2, mean)
b0 <- as.vector(f(cbind(y,x),0))
# Un autre échantillon de la même population
# A FAIRE : mettre ça dans une fonction
n <- 100
x0 <- rep(1,n)
x1 <- x0 + .1*rnorm(n)
x2 <- x0 + .1*rnorm(n)
x3 <- x0 + .1*rnorm(n)
y <- x0 + x1 + x2 + x3 + rnorm(n)
x <- cbind(x0,x1,x2,x3)
a1 <- x %*% b
a2 <- x %*% b0
boxplot( c(as.vector(a1), as.vector(a2)) ~
gl(2,n,2*n, c("ridge", "classic")),
horizontal=T )
title(main="Residues of a ridge regression")
A FAIRE
On peut utiliser le bootstrap pour valider un modèle de la manière suivante : on choisit le modèle, on le calcule sur des échantillons de bootstrap et on le valide sur l'échantillon initial. On a ainsi une estimation de l'erreur. Le problème, c'est qu'il y a des observations communes entre les échantillons de bootstrap et l'échantillon initial : cela va biaiser l'erreur vers 0.
On peut contourner le problème en calculant l'erreur sur l'observation numéro i uniquement à l'aide des échantillons de bootstrap qui ne la contiennent pas.
A FAIRE : en anglais, ça s'appelle "leave-out-one bootstrap", comment dit-on en français.
library(ipred)
Exemple d'utilisation du bootstrap : bande de confiance à 95% pour une régression avec des splines. Comparer avec les intervalles de prédiction fournis par la régression.
On vient de voir comment utiliser le bootstrap pour comparer des estimateurs, on peut aussi l'utiliser pour comparer des modèles, comme pour la validation croisée.
On peut rééchantillonner les couples (X,Y) On peut aussi rééchantillonner les résidus (normalisés comme il faut, de sorte que leur variance soit celle du bruit). plot( res ~ leverage )
A FAIRE : pour chaque exemple de bootstrap, faire des graphiques : boxplot, hist, qqnorm, jack.after.boot boot.ci to compute confidence intervals (several methods, some require very long computations).
A FAIRE : le maximum est un mauvais exemple -- regarder la distribution qu'on obtient. (et l'intervalle de confiance que je donne est stupide : il n'a aucune raison de s'étendre à gauche...)
A FAIRE : Classer tout ce qui suit...
Il y a un article sur le bootstrap et son implémentation dans R dans Rnews_2002-3.ps A FAIRE : URL library(help=bootstrap) Vieux library(help=boot) Plus récent Permet aussi de faire du bootstrap sur des données censurées à droite et des séries temporelles. ?censboot ?tsboot Also: multi-sample problems A FAIRE : influence function estimation from the bootstrap. ?empinf
Bootstrap et régression : souvent, on rééchantillone les résidus, pas les observations...
library(boot)
fit <- lm(calls ~ year, data=phones)
ph <- data.frame(phones, res=resid(fit), fitted=fitted(fit))
ph.fun <- function(data, i) {
d <- data
d$calls <- d$fitted + d$res[i]
coef(update(fit, data=d))
}
ph.lm.boot <- boot(ph, ph.fun, R=499)
ph.lm.boot
Les estimateurs que l'on obtient à l'aide de la méthode du maximum de vraissemblance (MLE) sont souvent biaisés. L'algorithme Jackknife, que voici, permet de transformer un tel estimateur U en un estimateur moins biaisé.
Diviser l'échantillon en N échantillons de même taille (par exemple, de taille 1) et calculer l'estimation Ui du paramètre. On pose alors
J = N U - (N-1)/N * (U1 + U2 + ... + Un)
Par exemple, si on lui donne la formule de la variance de la population (avec n, au dénominateur : c'est un estimateur biaisé) comme estimateur de la variance, il nous ressort la variance de l'échantillon (avec n-1 au dénominateur : il n'est plus biaisé).
A FAIRE : à vérifier
A FAIRE library(boot) ?jack.after.boot
Quand on ne connait pas du tout la loi d'une population, et si on a un échantillon, on peut utiliser cet échantillon pour faire des simulations concernant cette loi, tout simplement en faisant des tirages au hasard dans cet échantillon. C'est ce qu'on appelle la méthode de Monte-Carlo.
On peut l'utiliser pour comparer des estimateurs d'un paramètre de la loi, et ainsi choisir celui qui a la variance la plus faible.
A FAIRE library(boot) ?boot
A FAIRE :
utiliser le bootstrap pour : estimer la variance d'un estimateur calculer un intervalle de confiance faire des tests et calculer des p-valeurs
Exemple : prendre un échantillon, bootstrapper pour avoir un intervalle de confiance à 95% de la médiane. Représenter, en fonction de la taille de l'échantillon de bootstrap, la médiane médiane et les quantiles 0.025 et 0.975 de la médiane.
A FAIRE
A FAIRE :
# Utilisation du bootstrap (normal ou .632)ou de la validation # croisée pour estimer une erreur de classification. library(ipred) ?errortest xpdf ipred-examples.pdf
A FAIRE : cette partie est-elle bien à sa place ici ???
A FAIRE Graphe 1 : l'erreur de prédiction (sur l'échantillon
d'entrainement et sur l'échantillon de test) en fonction de la
complexité du modèle. Un modèle trop peu complexe ou trop complexe
donne de mauvaises prédictions : il faut savoir trouver le juste
milieu.
n <- 50
k <- 10
N <- 50
get.sample <- function (n) {
x <- runif(n)
y <- sin(2*pi*x)+.3*rnorm(n)
m <- data.frame(x=x, y=y)
m
}
error <- rep(NA,k)
for (i in 1:k) {
e <- rep(NA,N)
for (j in 1:N) {
training.sample <- get.sample(n)
x <- training.sample$x
y <- training.sample$y
r <- lm(y~poly(x,i))
test.sample <- get.sample(1000)
e[j] <- mean( (predict(r, data.frame(x=test.sample$x)) - test.sample$y)^2 )
}
error[i] <- mean(e)
}
plot(error, type='l', log='y',
xlab="Complexité du modèle",
ylab="Erreur")
On peut représenter le biais, la variance, l'erreur quadratique moyenne en fonction de la complexité du modèle, du paramètre de la ridge régression, du nombre de voisins dans la méthode des k plus proches voisins, de la vitesse de chute de la température pour le recuit simulé, du nombre de paramètres à retenir dans la régression sur les composantes principales ou les moindres carrés partiels, du degré dans une régression polynomiale, du nombre de noeuds dans une régression avec des splines cubiques restreintes, etc.
Il importe de ne pas confondre le biais, la variance et l'erreur quadratique moyenne.
A FAIRE : détailler
a <- 2
curve(dnorm(x-a), xlim=c(-3,5), ylim=c(0,.5))
abline(v=0, lty=3, col='red')
abline(v=a, lty=3, col='blue')
arrows(0,.45,a,.45, code=3, col='green')
arrows(a-pnorm(1), dnorm(pnorm(1)), a+pnorm(1), dnorm(pnorm(1)), code=3, col='red')
legend(par('usr')[2], par('usr')[4], xjust=1,
c("vraie valeur", "espérance de l'estimateur", "biais", "variance"),
lwd=1, lty=c(3,3,1,1),
col=c('red', 'blue', 'green', 'red') )
title(main="Biais et variance d'un estimateur biaisé")
Vincent Zoonekynd
<zoonek@math.jussieu.fr>
latest modification on Wed Oct 13 22:33:01 BST 2004