Sondage, et “prévision”

L’autre soir, @imparibus proposait sur son blog un billet passionnant sur l’élection présidentielle en France, avec des graphiques superbes, basé sur un lissage temporel des résultats des sondages pour le premier tour de l’élection présidentielle en France (qu’il a fait sous excel, on saluera la performance !). En reprenant les données du site http://www.lemonde.fr/, on peut obtenir la même chose assez simplement (j’ai repris – ou presque – les couleurs utilisées sur le site, François Hollande en rose, Nicolas Sarkozy en bleu, François Bayrou en orange, et Marine Le Pen en noir…). Pour commencer, le code en R pour importer les données et les manipuler est le suivant

sondage=read.table("http://freakonometrics.blog.free.fr/
public/data/sondage1ertour2012.csv",
header=TRUE,sep=";",dec=",")
sondage[36,1]="11/03/12"
sondeur=unique(sondage$SONDEUR)
sondage$date=as.Date(as.character(sondage$DATE),
"%d/%m/%y")

(je fais manuellement une correction de date car je me suis trompé en saisissant les chiffres à la main). A partir de là, on peut reprendre l’idée de faire une régression locale pour trouver une tendance,

CL=brewer.pal(6, "RdBu")
couleur=c(CL[1],CL[6],"orange","grey")
datefin=as.Date("2012/04/12")
k=1
plot(sondage$date,sondage[,k+2],col=couleur[k],
pch=19,cex=.7,xlab="",ylab="",ylim=c(0,40),
xlim=c(min(sondage$date),datefin))
rl=lowess(sondage$date,sondage[,k+2])
lines(rl,lwd=2,col=couleur[k])
for(k in 2:4){
points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7)
rl=lowess(sondage$date,sondage[,k+2])
lines(rl,lwd=2,col=couleur[k])
}

Le code a l’air long mais il faut définir les couleurs, générer une fenêtre graphique, mettre les régressions locales dedans, etc. En tant que telle, la commande qui permet de faire le lissage est juste

rl=lowess(sondage$date,sondage[,k+2])

Bon, on peut d’ailleurs faire toutes sortes de lissages, avec des splines par exemple (ce que j’ai davantage tendance à utiliser),

D=data.frame(date=min(sondage$date)+0:148)
library(splines)
k=1
plot(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prl=predict(rs,newdata=D)
prlse=sqrt(prl/100*(1-prl/100)/1000)*100
polygon(c(D$date,rev(D$date)),c(prl+2*prlse,
rev(prl-2*prlse)),col=CL[3],border=NA)
lines(D$date,prl,lwd=2,col=CL[1])
points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,)
for(k in 2:4){
points(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7)
}

avec, là encore, la commande suivante pour lisser

rs=lm(sondage[,k+2]~bs(date),data=sondage)

Cette fois le code est un peu plus long, parce que j’ai tracé un intervalle de confiance à 95% (intervalle de confiance classique, en supposant que 1000 personnes ont été interrogées, comme évoqué dans d’anciens billets),

http://freakonometrics.blog.free.fr/public/perso5/adic10.gif

On a ici la courbe suivante

avec l’intervalle de confiance pour François Hollande, mais on peut faire la même chose pour Nicolas Sarkozy,

k=2
plot(sondage$date,sondage[,k+2],
col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prl=predict(rs,newdata=D)
prlse=sqrt(prl/100*(1-prl/100)/1000)*100
polygon(c(D$date,rev(D$date)),c(prl+2*prlse,
rev(prl-2*prlse)),col=CL[4],border=NA)
lines(D$date,prl,lwd=2,col=CL[6])
points(sondage$date,sondage[,k+2],col=
couleur[k],pch=19,cex=.7,)
for(k in c(1,3:4)){
points(sondage$date,sondage[,k+2],col=
couleur[k],pch=19,cex=.7)
}

Mais comme auparavant, on peut aussi visualiser les “tendances” des quatre candidats en tête,

Voilà pour le début. En fait, idéalement, on voudrait faire un peu de prévision… Le hic est que les code usuel pour faire de la prévision (avec des ARIMA, du lissage exponentiel ou tout autre modèle classique) nécessitent de travailler avec des séries temporelles, telles qu’elles sont classiquement définies, c’est à dire “régulièrement espacées dans le temps“.

Je ne vais pas commencer à réclamer plus de sondages, loin de moins cette idée, mais pour mon modèle, j’avoue que j’aurais préféré avoir plus de points. Beaucoup plus de points. Un sondage par jour aurait été idéal en fait…

Qu’à cela ne tienne, on va simuler des sondages. L’idée est simple: on a une tendance qui nous donne une probabilité jour par jour (c’est exactement ce qui a été calculé pour tracer les courbes), via

D=data.frame(date=min(sondage$date)+0:148)
prl=predict(rs,newdata=D)

On va ensuite utiliser une hypothèse de normalité multivariée sur nos 4 candidats, auquel j’en ai en rajouté un cinquième fictif (mais ça ne servait à rien) en utilisant l’analogue multivariée de l’expression suivante (évoqué dans le premier billet de l’année sur les élections)

http://freakonometrics.blog.free.fr/public/perso5/ic-sondage01.gif

On peut alors générer, jour après jour, des sondages, en simulant des vecteurs Gaussiens. Je vais les générer indépendants les uns des autres, parce que c’est plus simple, et que cela ne me semble pas être une trop grosse hypothèse. Pour générer un ensemble de sondages, on utilise la fonction suivante

library(mnormt)
simulation=function(S){
proba=c(S,100-sum(S))/100
variance= -proba%*%t(proba)
diag(variance)=proba*(1-proba)
variance=variance/1000
return(rmnorm(1,proba[1:4],variance[1:4,1:4]))}
simulsondages=function(M){
Z=rep(NA,ncol(M))
for(i in 1:nrow(M)){
Z=rbind(Z,simulation(M[i,]))
}
return(Z[-1,])}
prediction4=matrix(NA,nrow(D),4)
for(k in 1:4){
rs=lm(sondage[,k+2]~bs(date),data=sondage)
prediction4[,k]=predict(rs,newdata=D)
}

Par exemple, si on la fait tourner une fois, on obtient le graphique suivant

set.seed(1)
S4=simulsondages(prediction4)
S100=100*S4
k=1
plot(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
}

(les courbes lissées sont celles obtenues sur les vrais sondages). Là on peut être heureux, parce qu’on a des vraies séries temporelles. On peut alors faire de la prévision, par exemple en faisant du lissage exponentiel (optimisé),

library(forecast)
k=1
X=S100[,k]
ETS=ets(X)
F=forecast(ETS,h=60)
fdate=max(D$date)+1:60
k=1
plot(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),lwd=2,col=couleur[k])
}
polygon(c(fdate,rev(fdate)),c(as.numeric(F$lower[,2]),
rev(as.numeric(F$upper[,2]))),col=CL[3],border=NA)
polygon(c(fdate,rev(fdate)),
c(as.numeric(F$lower[,1]),rev(as.numeric(F$upper[,1]))),
col=CL[2],border=NA)
lines(fdate,as.numeric(F$mean),lwd=2,col=CL[1])

Là encore le code peut paraître long, mais c’est surtout la partie associée au graphique qui prend de la place (par soucis purement esthétique, on passe du temps sur les codes graphiques depuis le début). Le code qui modélise la série, et qui la projette, est ici donné par les deux lignes suivantes

ETS=ets(X)
F=forecast(ETS,h=60)

Pour François Hollande, sur la simulation des sondages passés que l’on vient d’effectuer, on obtient la prévision suivante

On peut bien entendu faire une projection similaire pour Nicolas Sarkozy,

k=2
X=S100[,k]
ETS = ets(X)
F=forecast(ETS,h=60)
k=1
plot(D$date,S100[,k],col=couleur[k],pch=19,cex=.7,
xlab="",ylab="",ylim=c(0,40),xlim=
c(min(sondage$date),datefin))
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
for(k in 2:4){
points(D$date,S100[,k],col=couleur[k],
pch=19,cex=.7)
rs=lm(sondage[,k+2]~bs(date),data=sondage)
lines(D$date,predict(rs,newdata=D),
lwd=2,col=couleur[k])
}
polygon(c(fdate,rev(fdate)),c(as.numeric(F$lower[,2]),
rev(as.numeric(F$upper[,2]))),col=CL[4],border=NA)
polygon(c(fdate,rev(fdate)),
c(as.numeric(F$lower[,1]),rev(as.numeric(F$upper[,1]))),
col=CL[5],border=NA)
lines(fdate,as.numeric(F$mean),lwd=2,col=CL[6])

On notera que je ne me suis pas trop fatigué sur la projection: autant sur la génération de sondages passés on a tenu compte de corrélation (i.e. si un candidat a un score élevé, ça se fera au détriment des autres, d’où la corrélation négative utilisée dans la simulation), autant ici, les projections sont faites de manière complétement indépendantes. En se fatiguant un peu plus, on pourrait se lancer dans un modèle vectoriel gaussien. Ou mieux, on pourrait travailler sur des processus de Dirichlet (surtout que R a un package dédié à ce genre de modèle) parce qu’on travaille depuis le début sur des taux, comme évoqué auparavant. Mais commençons par faire simple pour un premier billet rapide sur ce sujet.

On peut ensuite s’amuser à générer plusieurs jeux de sondages, et de lancer des prévisions dessus,

set.seed(1)
for(sim in 1:20){
Ssim=simulsondages(prediction4)
F1=forecast(ets(100*Ssim[,1]),h=60)
F2=forecast(ets(100*Ssim[,2]),h=60)
lines(fdate,as.numeric(F1$mean),col=CL[1])
lines(fdate,as.numeric(F2$mean),col=CL[6])}

Une fois qu’on a fait tout ça, on a presque fini. On peut regarder au 22 avril qui est en tête (voire, par scénario, calculer la probabilité qu’un des deux candidats soit en tête), c’est à dire au soir du 1er tour

set.seed(1)
VICTOIRE=rep(NA,1000)
for(sim in 1:1000){
Ssim=simulsondages(prediction4)
F1=forecast(ets(100*Ssim[,1]),h=60)
F2=forecast(ets(100*Ssim[,2]),h=60)
VICTOIRE[sim]=(F1$mean[30]>F2$mean[30])}
mean(VICTOIRE)

Et voilà. Ah oui, je n’ai pas laissé le résultat. Tout d’abord parce que je suis un statisticien, pas unprédicateur. Mais aussi parce que si ce genre de pronostic amuse des gens… ils n’ont qu’à se mettre à R pour faire tourner les codes !

Publié dans Non classé | Marqué avec , | Laisser un commentaire

Le paradoxe de l’examen surprise…

La semaine dernière, Steve Landsburg reprenait (sur http://thebigquestions.com) sous une forme amusante un paradoxe que j’avais évoqué voilà deux ans sur le blog. Je ne résiste pas au plaisir de me replonger dans ce paradoxe…
Un professeur rédigeant son plan de cours veut mentionner un examen surprise, qui compterait pour 10% de la note finale.  Il annonce donc à ses étudiants, au premier cours, qu’un examen surprise (et un seul) se tiendrait à une date (surprise) parmi les 12 séances de cours. Mais comme le fait noter un élève, l’examen ne pourrait avoir lieu lors de la dernière séance (la 12ème), car cela signifierait qu’il n’y a pas eu d’examen les 11 premières séances, et donc l’examen surprise devrait forcément se tenir au 12ème cours (mais alors il ne serait plus surprise du tout). Bref, impossible d’avoir un examen surprise ce jour là. Mais dans ce cas, surenchérit une étudiante, il ne peut non plus avoir lieu à l’avant dernière séance (la 11ème), car là encore cela signifierait qu’il n’y a pas eu d’examen les 10 premières séances, et donc l’examen surprise devrait forcément au 11ème cours (car il ne peut pas avoir lieu au 12ème comme on vient de le voir, et là encore, il ne serait plus surprise du tout). Bref, de manière récurrente, on peut montrer que l’examen surprise a alors lieu forcément le premier jour. Mais là encore, il n’y a pas de surprise. Et donc impossible qu’il y ait un examen surprise si ce dernier a été annoncé. Moralité, pour faire un examen surprise, il ne faut surtout pas évoquer qu’il y en aura un. Mais comme on ne peut pas faire un examen non annoncé dans les modalités de contrôle des connaissances, il est impossible de faire des surprises aux étudiants à l’université… C’est triste, non ?

Publié dans Non classé | Marqué avec , , | Laisser un commentaire

Combien de temps profite-t-on de ses grands parents ?

Hier matin, je suis tombé un peu par hasard sur deux graphiques de l’INSEE (en France) avec l’age moyen des mères à l’accouchement, en fonction de rang de naissance de l’enfant, avec tout d’abord 1905-1965,

puis 1960-2000

Ces graphiques sont passionnants en soi – comme en ont témoigné pas mal de followers sur Twitter – mais ils m’ont fait m’interroger. En particulier, sur la croissance observée depuis 30 ans, qui me faisait penser à la tendance croissante observée sur les durées de vie. On n’a – malheureusement – pas accès au données complètes sur le site, mais on peut trouver d’autres données intéressantes (en l’occurrence l’age moyen à la naissance).

> agenaissance=read.table("http://freakonometrics.blog.free.fr/
public/data/agenaissance.csv",header=TRUE,sep=",")
> agenaissance$Age=as.character(agenaissance$AGE)
> agenaissance$AGE=as.numeric(substr(agenaissance$Age,1,2))+
+ as.numeric(substr(agenaissance$Age,4,4))/10
> plot(agenaissance$ANNEE+.5,agenaissance$AGE,
+ type="l",lwd=2,col="blue")

Visuellement, on retrouve la courbe en bleu foncée sur les graphiques ci-dessus,

On peut alors aller en cran plus loin, en se demandant non pas quel était l’âge moyen de la mère, mais de la grand-mère (au sens la mère de la mère)

> agenaissance$NAIS.MERE=(agenaissance$ANNEE+.5)-
+ agenaissance$AGE
> w=(trunc(agenaissance$NAIS.MERE-.5))
> rownames(agenaissance)=agenaissance$ANNEE
> a1=agenaissance[as.character(w),]$NAIS.MERE
> a2=agenaissance[as.character(w+1),]$NAIS.MERE
> p=agenaissance$NAIS.MERE-(w+.5)
> agenaissance$NAIS.GRD.MERE=(1-p)*a1+p*a2
> agenaissance$age.GRD.MERE=agenaissance$ANNEE+.5-
+ agenaissance$NAIS.GRD.MERE
> tail(agenaissance)
ANNEE  AGE   Age NAIS.MERE NAIS.GRD.MERE age.GRD.MERE
2000  2000 30.3 30,3     1970.2       1942.87        57.63
2001  2001 30.4 30,4     1971.1       1943.80        57.70
2002  2002 30.4 30,4     1972.1       1944.92        57.58
2003  2003 30.5 30,5     1973.0       1945.95        57.55
2004  2004 30.5 30,5     1974.0       1947.05        57.45
2005  2005 30.6 30,6     1974.9       1948.04        57.46
> plot(agenaissance$ANNEE+.5,agenaissance$age.GRD.MERE,
+ type="l",lwd=2,col="red")

Là encore, on peut visualiser l’âge de la grand-mère maternelle à la naissance

A partir de là, on peut se demander combien de temps on profite de ses grands-parents (ou tout du moins ici de sa grand mère maternelle), en se basant sur les calculs d’espérance de vie résiduelle. En utilisant le modèle de Lee-Carter pour modéliser les taux de décès annuel, et en extrapolant sur le siècle en cours, on peut extrapoler les espérances de vie résiduelles.

> Deces <- read.table("http://freakonometrics.free.fr/
Deces-France.txt",header=TRUE)
> Expo  <- read.table("http://freakonometrics.free.fr/
Exposures-France.txt",header=TRUE,skip=2)
> Deces$Age <- as.numeric(as.character(Deces$Age))
> Deces$Age[is.na(Deces$Age)] <- 110
> Expo$Age <- as.numeric(as.character(Expo$Age))
> Expo$Age[is.na(Expo$Age)] <- 110
>  library(forecast)
>  library(demography)
>  YEAR <- unique(Deces$Year);nC=length(YEAR)
>  AGE  <- unique(Deces$Age);nL=length(AGE)
>  MUF  <- matrix(Deces$Female/Expo$Female,nL,nC)
>  POPF <- matrix(Expo$Female,nL,nC)
>  BASEF <- demogdata(data=MUF, pop=POPF,ages=AGE,
+ years=YEAR, type="mortality",
+  label="France", name="Femmes", lambda=1)
> LCF <- lca(BASEF)
> LCFf<-forecast(LCF,h=100)
> A <- LCF$ax
> B <- LCF$bx
> K1 <- LCF$kt
> K2 <- K1[length(K1)]+LCFf$kt.f$mean
> K <- c(K1,K2)
> MU <- matrix(NA,length(A),length(K))
> for(i in 1:length(A)){
+ for(j in 1:length(K)){
+ MU[i,j] <- exp(A[i]+B[i]*K[j]) }}
> esp.vie = function(xentier,T){
+ s <- seq(0,99-xentier-1)
+ MUd <- MU[xentier+1+s,T+s-1898]
+ Pxt <- cumprod(exp(-diag(MUd)))
+ ext <- sum(Pxt)
+ return(ext) }
> EVIE = function(x,T){
+ x1 <- trunc(x)
+ x2 <- x1+1
+ return((1-(x-x1))*esp.vie(x1,T)+(x-x1)*esp.vie(x2,T)) }
> agenaissance$EV=NA
> for(i in 1:100){
+ t <- 2006-i
+ agenaissance$EV[agenaissance$ANNEE==t]=
+ EVIE(x=agenaissance$age.GRD.MERE[
+ agenaissance$ANNEE==t],t) }
> tail(agenaissance)
ANNEE  AGE   Age NAIS.MERE NAIS.GRD.MERE age.GRD.MERE       EV
2000 30.3 30,3     1970.2       1942.87        57.63 29.13876
2001 30.4 30,4     1971.1       1943.80        57.70 29.17047
2002 30.4 30,4     1972.1       1944.92        57.58 29.39027
2003 30.5 30,5     1973.0       1945.95        57.55 29.52041
2004 30.5 30,5     1974.0       1947.05        57.45 29.72511
2005 30.6 30,6     1974.9       1948.04        57.46 29.80398

Autrement dit, sur la dernière ligne, l’espérance de vie (résiduelle) pour une femme de 57.46 ans en 2005 était d’environ 29.80 ans. On peut alors visualiser non seulement l’âge moyen de sa grand-mère à la naissance, mais son espérance de vie résiduelle,

> plot(agenaissance$ANNEE+.5,agenaissance$EV,
+ type="l",lwd=2,col="purple")

On note que depuis 30 ans, en France, la durée (moyenne) pendant laquelle les petits-enfants vont profiter de leur grands parents s’est stabilisé à une trentaine d’années. On peut aussi continuer, et remonter d’un cran (en refaisant tourner le code avec quelques modifications): on a alors l’âge (moyen) de son arrière grand mère à la naissance,

et la durée de vie (résiduelle) des arrière grand mères

On manque ici de données, mais il semble que l’on profite – en moyenne – environ 5 ans de son arrière grand mère. Maintenant on peut aussi s’interroger sur les limites de cette étude rapide. En particulier, de même qu’il existe une corrélation forte entre les durées de vie de conjoints (e.g. broken heart syndrom de Jagger & Sutton (1991)), on peut se demander si la naissance d’enfants et de petits-enfants a un impact sur la durée de vie résiduelle d’une personne (ou si on peut supposer l’indépendance comme on l’a fait ici).

Publié dans Non classé | Laisser un commentaire

Plus intéressants que les nobels: les ig-nobels

Il y a quelques semaines, Alessandro Pluchino, Andrea Rapisarda et Cesare Garofalo ont obtenu un des IG nobels (ici, probablement la plus haute distinction scientifique qui soit) pour avoir proposé une étude numérique du principe de Peter. Comme le notait Rue89, ils ont fait la “démonstration mathématique que la promotion aléatoire permet de gagner en efficacité dans une entreprise” ().

  • Le principe de Peter

Historiquement, tout vient d’un livre de Laurence Peter et Raymond Hull, paru en 1969, et intitulé précisément Peter principle. On y retrouve une formulation connu en management, avec le temps, tout poste sera occupé par un incompétent incapable d’en assumer la responsabilité“. Pour cela, on part de principes très simples, et que tout le monde admettra. Le premier est que quelqu’un de compétent à un poste est promu (à un niveau supérieur). Le second est que quelqu’un d’incompétent à un poste n’est pas promu à un niveau supérieur.

Aussi, un employé ne peut pas rester à un poste où il est compétent (puisqu’il est promu) et à l  suite des promotions, il devrait finir par atteindre un poste auquel il sera incompétent, il ne recevra plus de promotion, il restera donc indéfiniment à ce poste (pour lequel il est incompétent). Moralité, tous les postes finissent par être occupés par des employés incompétents, et à un instant donné on arrive à la formule lapidaire annonçant que la majorité du travail est effectuée par des salariés n’ayant pas encore atteint leur “seuil d’incompétence“.

Bon, c’est un peu simpliste, supposant qu’il existe une relation d’ordre, en particulier qu’une personne bonne à un niveau donné serait forcément bonne aux niveaux “inférieurs”. Il existe en particulier un compétence définie en une dimension unique, et donc avec un ordre naturel,

  • Le papier de Pluchino, Rapisarda et Garofalo

Le papier d’Alessandro Pluchino, Andrea Rapisarda et Cesare Garofalo (en ligne ici) essayent de tester numériquement ce principe, et le valident. Ils font en particulier un test très intéressant sur la stratégie des promotions. Ils montrent – sous quelques hypothèses de modélisation – que la stratégie visant à promouvoir le “meilleur” n’ait pas forcément optimale, et qu’on peut avoir intérêt à avoir une stratégie mixte en retenant parfois le plus mauvais… avec un seuil de mélange proche de 50%. Cette stratégie peut aussi se voir comme un stratégie simple: soit on offre une promotion au meilleur,  soit on offre une promotion par tirage au sort…. et comme c’est un peu ce qui se passe dans les tournois de soccer (si on pense que s’en remettre au tir au but c’est s’en remettre au hasard) on a un moyen simple pour tester cette théorie….

  • Étude sur les matchs de coupe du monde de soccer

J’ai regardé les résultats sur les coupes du monde allant de 1982 à 2010, en phase finale, soit 8+4+2=14 matchs par coupe, sauf en 1982, où je n’ai que deux matchs. Soit 100 matchs en tout. A chaque match, soit le gagnant s’est fait sur le score (le vainqueur gagne) soit sur le hasard (aux tirs au buts). Et pour tous ces matchs, je regarde si l’équipe qualifié a gagné, ou pas, le match suivant. J’obtiens la matrice suivante,

matchs gagnant perdant total
au mérite 46 34 80
au hasard 4 16 20
total 50 50 100

Si on fait un test d’indépendance du chi-deux sur ce tableau de contingence, on obtient qu’il n’y pas indépendance entre les variables,
> M=matrix(c(46,4,34,16),2,2)
> chisq.test(M)
Pearson’s Chi-squared test with Yates’ continuity correction
data:  M
X-squared = 7.5625, df = 1, p-value = 0.00596
Autrement dit, le fait d’être gagnant – ou perdant – à un match (entre les quarts et la finale) dépend fortement de la manière dont on s’est qualifié… Promouvoir en faisant un tirage au sort ne semble donc pas vraiment optimal, l’équipe qualifiée par hasard ne passant souvent pas l’étape suivante… Mais je suis intéressé par toute analyse sportive de ce principe, donc si quelqu’un a des pistes, je suis preneur ! (car les données ne manquent pas)

Publié dans Non classé | Marqué avec , | Laisser un commentaire

70% des salariés français touchent moins de 1 500 € par mois ?

Hier, @ me demandait de commenter l’information comme quoi “70% des salariés français touchaient moins de 1 500 € par mois“. Tout d’abord, rappelons que selon l’INSEE, en 2008, le salaire médian net en France est de l’ordre de 1655 € par mois pour un temps plein (cf ici par exemple). Comme c’est la médiane (je renvoie ici pour une bonne blague du nouvel obs sur la médiane, comme quoi il ne faut pas toujours croire ce qu’on lit), c’est que 50% des français touchent ce salaire, ou plutôt toucheraient ce salaire s’il travaillaient à temps plein…. Le hic est que généralement l’insee normalise les données, et ne parle que de taux plein…

Mais heureusement l’insee met aussi a disposition des bases, en particulier la base dads2007, ici, avec un fichier de 25Mo (zippé) , contenant près de 2 millions d’observations, dont le temps travaillé en 2007 et le salaire (ou plutôt une tranche de salaire). Le code est plutôt simple (même si j’ai travaillé sur la première moitié de la base, suite à des soucis d’importation sur mon petit portable),

> library(foreign)
> base=read.dbf("D:\\r-data\\salaries07.dbf")
> TRANCHE=c(0,100,300,600,1000,1300,2000,3000,
+ 4000,5000,7000,9000,11000,12500,14000,16000,
+ 18000,20000,24000,30000,1000000000)
> salinf=TRANCHE[1:20]
> salsup=TRANCHE[2:21]

Bref,
on peut récupérer des salaires médians, par temps plein (afin de voir si nos données sont représentatives), en notant qu’un temps plein c’est 1820 heures dans l’année (cf ici)

> Xinf=salinf[base$TRNNETO+1]/base$NBHEUR_TOT*1820/12
> Xsup=salsup[base$TRNNETO+1]/base$NBHEUR_TOT*1820/12
> quantile(Xinf,.5,na.rm=TRUE)
 50%
1484.745
> quantile(Xsup,.5,na.rm=TRUE)
 50%
1685.185

Comme nous ne sommes pas trop loin des 1655 indiqués, j’aurais tendance à faire confiance à nos chiffres (en se plaçant proches de la borne supérieure d’ailleurs). Mais si on travaille sur le revenu, et pas le revenu pour un temps plein, les chiffres sont un peu différents,

> Xinf=salinf[base$TRNNETO+1]/12
> Xsup=salsup[base$TRNNETO+1]/12
> quantile(Xinf,.5,na.rm=TRUE)
 50%
1333.333
> quantile(Xsup,.5,na.rm=TRUE)
 50%
1500

Pour revenir à nos moutons, si on s’intéresse à la proportion des français qui touchent moins de 1500 € par mois, j’obtiens

> mean(Xinf<=1500)
[1] 0.6118895
> mean(Xsup<=1500)
[1] 0.5236287

C’est à dire qu’à la louche, je pense que 55% des français touchent moins de 1500 € par mois. A condition que les chiffres soient en net. Si la citation parlait du brut, ça correspond en gros à du 1200 € par mois, et on a alors les proportions suivantes,

> mean(Xinf<=1500*.8)
[1] 0.4215607
> mean(Xsup<=1500*.8)
[1] 0.3143752

D’ailleurs un petite remarque sur ces histoires de revenu et de temps de travail. Si le temps travaillé par mois et le salaire horaire étaient des variables indépendantes, j’aurais pu faire des calculs rapides à partir de la distribution du temps de travail. Malheureusement, ces deux données ne sont pas indépendantes. Pour s’en convaincre, calculons le produit des espérances (salaire moyen horaire par temps moyen travaillé par mois) et l’espérance du produit (revenu moyen)

> Xinf=salinf[base$TRNNETO+1]/base$NBHEUR_TOT*1820/12
> Xsup=salsup[base$TRNNETO+1]/base$NBHEUR_TOT*1820/12
> mean(Xinf[Xinf<1000000000000],na.rm=TRUE)*
+ mean(base$NBHEUR_TOT/1820)
[1] 1609.669
> mean(salinf[base$TRNNETO+1])
[1] 16832.14

qui ne sont manifestement pas égaux. On a même un dépendance positive entre le salaire horaire et le temps de travail mensuel.  On retrouve ici le fait que les personnes qui effectuent les travaux avec les salaires horaires les plus bas sont aussi celles qui travaillent le moins. D’un point de vue technique, pour étudier les inégalités temps de travail/revenu horaire, c’est un peu compliqué à cause du produit. C’est pour ça qu’avec Stéphane Mussard on avait introduit la décomposition  des mesures d’inégalités pour des inégalités multiplicatives, ici. Mais on sort du cadre du papier, et ma réponse serait que non, je doute que la personne soit de bonne foi (la seule source que j’ai trouvé sur internet donnant ces chiffres était ici).

Publié dans Non classé | Marqué avec | Laisser un commentaire

MacGyver et la théorie des groupes

L’autre jour, quelqu’un me félicitait pour mon blog, c’était une étudiante en maths, qui disait comprendre certaines choses en probas et en stats (entre autres sur les applications possibles)…. mais pour la théorie des groupes, elle ne voyait toujours pas d’intérêt…
J’étais bien entendu ravi que mon blog serve à quelque chose, mais en fait la théorie des groupes peut avoir des applications amusantes (bon, tout dépend aussi de ses centres d’amusement). Par exemple MacGyver pourrait avoir eu besoin de connaître des résultats sur les groupes
pour sortir des prisonniers des geôles afghanes…

Considérons le gardien de prison sadique qui annonce le “jeu” suivant (Anatole, un collègue ici me l’avait suggéré l’autre jour en me mettant au défi de trouver une solution…. mais faut pas pousser, je suis plus fort que MacGyver ! bon, et il faut dire que le jeu est classique…on le trouve expliqué ici ou ):
Tous les prisonniers (disons qu’ils sont 20) se voient attribuer un numéro,
entre 1 et 20. Le premier prisonnier est appelé et se retrouve face à 20 portes, derrière lesquelles sont cachés des numéros (allant de 1 à 20 là encore). Il a droit d’en ouvrir 10 (au maximum) et il doit retrouver son numéro. S’il n’y arrive pas,tous les prisonniers sont exécutés (lui compris). S’il y arrive, il n’est pas sauvé pour autant car le commandant appelle le second prisonnier, qui se retrouve à son tour devant les 20 portes, et doit trouver son numéro… Etc.

Les prisonniers commencent à paniquer, car tout le monde doit trouver son numéro pour être sauvé. En raisonnant rapidement, un des prisonniers explique qu’un seul prisonnier aurait 50% de chances d’y arriver. Avec deux, comme ils ont tous les deux 50% de chances, on
tombe à 25%, etc. Avec n prisonniers, on est à 1/2n, qui est de l’ordre d’une chance sur un million avec 20 prisonniers…. Et qu’ils ont de la chance de n’être que 20 ! Avec 30 prisonniers, ils
auraient une chance sur un milliard de s’en sortir (ce qui est faible, on en conviendra).
Mais MacGyver, qui avait pris le cours Théorie des Groupes 101 à l’université (et qui avait R installé dans sa montre digitale) dit

  • pas du tout, j’ai une stratégie qui nous donne environ une chance sur trois
    de tous s’en sortir (vu comme ça c’est pas terrible, mais par rapport à
    1 chance sur un milliard, c’est vachement bien),
  • et accessoirement, cette probabilité ne dépend pas du nombre de prisonnier (tant que l’on puisse ouvrir moitié moins de portes qu’il y a de prisonniers)

Bon, sur le moment, les autres prisonniers se fichent du fait que même s’ils étaient 1000 ils auraient une chance sur trois de survivre, ils veulent comprendre la stratégie de Mac… Comme le rappelle Mac, le commandant ne fait que choisir une permutation de http://perso.univ-rennes1.fr/arthur.charpentier/latex/prison99.png, même s’il est sadique (ou parce qu’il est sadique, les profs de maths sont sadiques, tout le monde sait ça). Et toute permutation peut se décomposer en un produit ds combinaisons cycliques (on parle de la décomposition en produit de cycles à supports disjoints, c’est un truc utilise je crois pour les amateurs de Rubik’s cube),
Par exemple, considérons la permutation suivante (c’est la première que j’ai eu sous R),

http://perso.univ-rennes1.fr/arthur.charpentier/latex/prison01.png

Elle se décompose en 5 permutations cycliques

http://perso.univ-rennes1.fr/arthur.charpentier/latex/prison02.png

La plus grande est la dernière, et elle permute sur 7 nombres. Autrement dit, si chacun utilise cette stratégie, il ouvrira au plus 7 portes. Le point important est qu’ils ne peuvent pas être pris dans
une autre permutation que la leur, car ils ont pris leur porte, donc au pire, il tire le numéro de la première porte qu’ils ont ouvert… qui est le numéro qu’ils doivent tirer pour survivre !
Regardons sur notre exemple
> set.seed(1)
> portes = sample(1:n,size=n,replace=FALSE)
>  prisonnier=1
>  essai = c(prisonnier)
>  nessai = 1
>  i = porte[prisonnier]
>  while(nessai <= 50) {
+  essai = c(essai, i)
+  if(i == prisonnier) {
+  cat(“Le prisonnier”, prisonnier, “a ouvert”, nessai, “portes\n”)
+  break; }
+  else {
+  i = porte[i]}
+ n essai =n essai+1 }
Le prisonnier 1 a ouvert 4 portes
>  cat(“Le prisonnier”, prisonnier, “a ouvert les portes suivantes”, paste(essai[-length(essai)], collapse=” – “), “\n”)
Le prisonnier 1 a ouvert les portes suivantes 1 – 6 – 14 – 10

soit ici
Le prisonnier 1 a ouvert 4 portes
Le prisonnier 1 a ouvert les portes suivantes 1 – 6 – 14 – 10
Le prisonnier 2 a ouvert 7 portes
Le prisonnier 2 a ouvert les portes suivantes 2 – 8 – 9 – 19 – 18 – 17 – 12
Le prisonnier 3 a ouvert 2 portes
Le prisonnier 3 a ouvert les portes suivantes 3 – 11
Le prisonnier 4 a ouvert 5 portes
Le prisonnier 4 a ouvert les portes suivantes 4 – 16 – 7 – 15 – 5
Le prisonnier 5 a ouvert 5 portes
Le prisonnier 5 a ouvert les portes suivantes 5 – 4 – 16 – 7 – 15
Le prisonnier 6 a ouvert 4 portes
Le prisonnier 6 a ouvert les portes suivantes 6 – 14 – 10 – 1
Le prisonnier 7 a ouvert 5 portes
Le prisonnier 7 a ouvert les portes suivantes 7 – 15 – 5 – 4 – 16
Le prisonnier 8 a ouvert 7 portes
Le prisonnier 8 a ouvert les portes suivantes 8 – 9 – 19 – 18 – 17 – 12 – 2
Le prisonnier 9 a ouvert 7 portes
Le prisonnier 9 a ouvert les portes suivantes 9 – 19 – 18 – 17 – 12 – 2 – 8
Le prisonnier 10 a ouvert 4 portes
Le prisonnier 10 a ouvert les portes suivantes 10 – 1 – 6 – 14
Le prisonnier 11 a ouvert 2 portes
Le prisonnier 11 a ouvert les portes suivantes 11 – 3
Le prisonnier 12 a ouvert 7 portes
Le prisonnier 12 a ouvert les portes suivantes 12 – 2 – 8 – 9 – 19 – 18 – 17
Le prisonnier 13 a ouvert 2 portes
Le prisonnier 13 a ouvert les portes suivantes 13 – 20
Le prisonnier 14 a ouvert 4 portes
Le prisonnier 14 a ouvert les portes suivantes 14 – 10 – 1 – 6
Le prisonnier 15 a ouvert 5 portes
Le prisonnier 15 a ouvert les portes suivantes 15 – 5 – 4 – 16 – 7
Le prisonnier 16 a ouvert 5 portes
Le prisonnier 16 a ouvert les portes suivantes 16 – 7 – 15 – 5 – 4
Le prisonnier 17 a ouvert 7 portes
Le prisonnier 17 a ouvert les portes suivantes 17 – 12 – 2 – 8 – 9 – 19 – 18
Le prisonnier 18 a ouvert 7 portes
Le prisonnier 18 a ouvert les portes suivantes 18 – 17 – 12 – 2 – 8 – 9 – 19
Le prisonnier 19 a ouvert 7 portes
Le prisonnier 19 a ouvert les portes suivantes 19 – 18 – 17 – 12 – 2 – 8 – 9
Le prisonnier 20 a ouvert 2 portes
Le prisonnier 20 a ouvert les portes suivantes 20 – 13
En adaptant un code développé ici par Matt Asher, on peut d’ailleurs visualiser tout ça sur un petit dessin,
,

Bref, la grande question est “quelle est la probabilité que la plus grande sous permutation cyclique soit sur un ensemble de taille supérieure à 10
? On peut faire un peu de Monte Carlo pour voir…. c’est ce que fait Matt. Mais sinon on peut aussi faire du dénombrement (ce qui est pareil que les probas). Le nombre de permutations contenant un cycle de longueur supérieure à 10 (strictement) se calcule de la manière suivante: soit k la taille du cycle, alors je peux choisir les http://perso.univ-rennes1.fr/arthur.charpentier/latex/prison00.png nombres de http://perso.univ-rennes1.fr/arthur.charpentier/latex/prison03.png manières. Il y a alors http://perso.univ-rennes1.fr/arthur.charpentier/latex/prison04.png permutations circulaires de ces k nombres, et on permute les http://perso.univ-rennes1.fr/arthur.charpentier/latex/prison05.png nombres qui reste comme on veut, ce qui donne http://perso.univ-rennes1.fr/arthur.charpentier/latex/prison06.png possibilités. Bref, j’en suis à

http://perso.univ-rennes1.fr/arthur.charpentier/latex/prison07.png

permutations ayant un cycle de longueur k, soit en tout,

http://perso.univ-rennes1.fr/arthur.charpentier/latex/prison0_.png

Or comme on a en tout http://perso.univ-rennes1.fr/arthur.charpentier/latex/prison08.png permutations possibles, la probabilité de s’en sortir vivant est

http://perso.univ-rennes1.fr/arthur.charpentier/latex/prison09.png

Et là où Mac est fort, c’est qu’effectivement, avec n prisonniers

http://perso.univ-rennes1.fr/arthur.charpentier/latex/prison10.png

Soit environ 30% de chances de s’en sortir….
Moralité, la théorie des groupes peut être marrante et intéressante…. Bon, c’est la seule application à la vraie vie qui me vient à l’esprit, et à moins d’être enfermé en prison avec un geôlier sadique, c’est dur de trouver des applications. N’empêche que la solution est passionnante (mais MacGyver a toujours été mon héros pour ça…)

Publié dans Non classé | Marqué avec , , | Laisser un commentaire

Des étés pluvieux en Bretagne ? une réalité statistique…

Pour compléter le précédant billet (ici) on peut se demander en quoi la Bretagne est différente des autres régions françaises… Nous avions vu ici le niveau de précipitation moyen, jour après jours pendant les mois d’été, en Bretagne. A Rennes. En revanche, à Paris on obtient la moyenne suivante,

que l’on peut comparer à Marseille,

ou encore à Strasbourg,

Sur la figure ci-dessous, on voit que la probabilité d’avoir de la pluie à Paris (au sens au moins 0.1 mm d’eau dans la journée, en trait gras bleu, au moins 2 mm d’eau dans la journée, en trait bleu) est supérieure à la probabilité d’avoir de la pluie à Rennes (respectivement en bleu clair gras, et en bleu clair fin)

On est certes très au dessus de Marseille,

mais très en dessous de Strasbourg,

Mais au delà des lois marginales, ces villes sont différentes de la Bretagne si l’on regarde les matrices de transition.

  • Transition d’un jour sur l’autre

Pour Rennes, si on regarde jour après jour, on obtient

jour 
beau
temps
pluie
jour 
beau temps
1955 612 2567
pluie
606 723 1329

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
76,15 % 23,85 %
pluie
45,60 % 54,40 %

Pour Paris, la probabilité de transition jour après jour a la forme suivante

jour 
beau
temps
pluie
jour 
beau temps
2689 959 3648
pluie
946 1466 2412

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
73,71 % 26,29 %
pluie
39,22 % 60,78 %

Pour Marseille, la probabilité de transition jour après jour a la forme suivante

jour 
beau
temps
pluie
jour 
beau temps
2527 375 2902
pluie
362 216 578

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
87,08 % 12,92 %
pluie
62,63 % 37,37 %

Pour Strasbourg, la probabilité de transition jour après jour a la forme suivante

jour 
beau
temps
pluie
jour 
beau temps
31 128 159
pluie
132 1464 1596

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
19,50 % 80,50 %
pluie
8,27 % 91,73 %
  • Transition d’une semaine sur l’autre

Si en revanche on regarde les matrices de transition semaine par semaine, on a des résultats assez différents. Une bonne semaine signifie aucun jour avec plus de 2 mm de pluie.
Pour Rennes, si on regarde semaine après semaine

jour 
beau
temps
pluie
jour 
beau temps
379 25 404
pluie
26 7 33

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
93,81 % 6,19 %
pluie
78,79 % 21,21 %

Pour Paris, la probabilité de transition semaine après semaine a la forme suivante

jour 
beau
temps
pluie
jour 
beau temps
576 46 622
pluie
53 4 57

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
92,60 % 7,40 %
pluie
92,98 % 7,02 %

Pour Marseille, la probabilité de transition semaine après semaine a la forme suivante

jour 
beau
temps
pluie
jour 
beau temps
274 59 333
pluie
47 9 56

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
82,28 % 17,72 %
pluie
83,93 % 16,07 %

Pour Strasbourg, la probabilité de transition semaine après semaine a la forme suivante

jour 
beau
temps
pluie
jour 
beau temps
1494 614 2018
pluie
613 939 1552

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
70,87 % 29,13 %
pluie
39,50 % 60,50 %

Les tests du chi deux, d’indépendance d’une semaine sur l’autre donnent

  • à Rennes, une statistique du chi-deux de 8,054, soit une p-value de 0,45%
  • à Paris, une statistique du chi-deux de 0,025, soit une p-value de 87,26%
  • à Marseille, une statistique du chi-deux de 0,012, soit une p-value de 91,24%
  • à Strasbourg, une statistique du chi-deux de 0,7649, soit une p-value de 38,18%

autrement dit l’hypothèse d’indépendance est acceptée partout, sauf à Rennes….

  • Moralité ?

De manière assez paradoxale, on prétend que la Bretagne a un temps changeant, et pour reprendre le titre du précédant billet, effectivement, en Bretagne, il peut faire beau plusieurs fois par jour. Mais sur le long terme, d’une semaine sur l’autre, le temps est au contraire très corrélé, contrairement aux autres régions. A Paris, Marseille ou Strasbourg, qu’il ait fait beau ou qu’il ait plus la semaine précédente, cela n’apporte aucune information sur la probabilité d’avoir de la pluie la semaine où l’on vient en vacances…. Mais pas en Bretagne: manifestement, il existe donc des étés pourris, où il pourra pleuvoir toutes les semaines, et des étés superbes où il ne pleut jamais….

Publié dans Non classé | Marqué avec , , , , , , , | Laisser un commentaire

En Bretagne, il fait beau plusieurs fois par jour

bon, et à l’occasion il peut pleuvoir un peu…. Il peut donc être intéressant pour planifier un peu ses vacances de calculer la probabilité d’avoir de la pluie.

Les calculs ont été fait sur les données de pluviométrie à Rennes, en ligne sur le site eca&d, ici (données de qualité gratuites).

  • Probabilité d’avoir de la pluie, pendant les vacances d’été

Considérons ici la série du niveau de précipitation par jour, en par 0,1 mm,

et la probabilité d’avoir de la pluie dans la journée (régression logistique),

Bref, il peut pleuvoir à Rennes l’été. Mais ça n’aide pas vraiment pour planifier ses vacances. Car si ça se trouve, il y a des étés sans pluie, et des étés où il pleut sans cesse.

  • Dynamique et matrice de transition (par jour et par semaine)

Autrement dit, au lieu de regarder les lois marginales (comme la probabilité d’avoir de la pluie dans la journée), on peut s’intéresser à la dynamique de la série, modélisée sous la forme d’une chaîne de Markov. Si on regarde jour après jour, avec les 30 mois de juillet et août entre 1980 et 2009,

jour 
beau
temps
pluie
jour 
beau temps
871 298 1169
pluie
292 335 627

ce qui donne les probabilités de transition suivantes,

jour 
beau
temps
pluie
jour 
beau temps
75,51 % 25,49 %
pluie
46,57 % 53, 43 %

autrement, s’il fait beau aujourd’hui, on a 3 chances sur 4 d’avoir du beau temps demain.
Si on regarde semaine après semaine, où l’intérêt sont les semaines sans pluie

,

semaine 
beau
temps
pluie
semaine 
beau temps
13 26 39
pluie
23 140 163

avec là aussi les probabilités de transition suivantes,

semaine 
beau
temps
pluie
semaine 
beau temps
33,33 % 66,67 %
pluie
14,11 % 85,89 %

Si on regarde semaine après semaine, où l’intérêt sont les semaines avec six jours sans pluie (on s’autorise une journée de pluie),

semaine 
beau
temps
pluie
semaine 
beau temps
36 36 72
pluie
30 100 130

ce qui donne les probabilités de transition suivantes,

semaine 
beau
temps
pluie
semaine 
beau temps
50,00 % 50,00 %
pluie
23,08 % 76,92 %

Bon, maintenant, si on pense qu’avoir 2 mm de pluie dans la journée, c’est juste un peu d’humidité dans l’air, les matrices de transitions sont sensiblement différentes, dans le cas où on s’autorise une journée dans la semaine avec un peu d’humidité pour parler de beau temps,

semaine 
beau
temps
pluie
semaine 
beau temps
106 36 142
pluie
31 29 60

ce qui donne les probabilités de transition suivantes,

semaine 
beau
temps
pluie
semaine 
beau temps
75,65 % 25,35 %
pluie
51,67 % 48,33 %

Autrement dit, on retrouve une matrice proche de celle que nous avions sur les données journalière: s’il a fait beau la semaine avant de venir en Bretagne, on a 3 chances sur 4 d’avoir du beau temps. En revanche, s’il a fait mauvais, on a une chance sur deux d’avoir beau temps la semaine suivante.
Si notre définition de pluie est encore plus laxiste (il faut qu’il y ait eu un déluge une journée dans la semaine, à savoir plus de 2 cm d’eau dans la journée), alors cette fois, on obtient,

semaine 
beau
temps
pluie
semaine 
beau temps
171 12 183
pluie
14 5 130

ce qui donne les probabilités de transition suivantes,

semaine 
beau
temps
pluie
semaine 
beau temps
93,44 % 6,56 %
pluie
73,68 % 26,32 %

Autrement dit, s’il y a eu au moins une très mauvaise journée la semaine passée, on a 1 chance sur 4 d’en avoir également une la semaine suivante. En revanche, s’il a fait beau tout le temps, on a presque 95% de chances d’avoir du beau temps la semaine suivante. Moralité, quelle que soit la définition retenue pour définir le beau temps, le temps à Rennes n’est pas indépendant d’une semaine sur l’autre: il y a manifestement des étés pluvieux, et d’autre non.

Publié dans Non classé | Marqué avec , , , , , | Laisser un commentaire