Tag Archives: statistique

De la symbolique des grandeurs statistiques

Ce matin, je suis suis un peu emporté sur Twitter, alors que d’ordinaire je préfère éviter…

Il faut dire que j’étais tombé sur un article dans les Echos qui m’a énervé…

J’ai ensuite passé beaucoup de temps à essayer de comprendre pourquoi j’avais été autant énervé par ce titre, cet article… car il n’apprenait rien de bien nouveau… Je pense que c’est vraiment la forme qui m’a agacé, et la symbolique du” temps” y joue (je pense) pour beaucoup. Pourtant, utiliser l’année calendaire comme métaphore pour expliquer des statistiques n’est pas nouveau.

Je pense que la première fois que je ai vu cette idée, je devais avoir dix ans, quand on nous parlait de l’histoire de l’univers, de la terre et l’humanité. Parce que mine de rien, manipuler des grands chiffres n’aide pas à avoir une représentation (ce qui est le cas quand on parle du budget de l’Etat). Oui, le big-bang a eu lieu il y a un peu moins de 15 milliards d’années, Lucy date d’environ 3,2 millions d’années, et le Christ est né il y a 2000 ans. Mais jongler entre les milliers, les millions et les milions est compliqué. Alors dans les cours d’histoire, on utilise une image basée sur le calendrier. Le Big Bang a eu lieu le 1er janvier à zéro heure et nous sommes aujourd’hui le 31 décembre à minuit, autrement dit chaque seconde équivaut à un peu plus de 400 années. Début septembre la terre se forme, les poissons et les premières plantes terrestres arrivent juste avant l’hiver, vers le 20 décembre. Lucy se met debout à 22h30, les Pyramides d’Égypte datent de minuit moins dix secondes, et le Christ naît à minuit moins cinq secondes. Voilà en gros pour l’image. En 90 minutes (sur une année) on a l’histoire de l’humanité, et en 10 secondes, l’histoire telle qu’on l’apprend à l’école. Mais dans cette image, on transforme du temps en un autre temps. Un peu comme dans la métaphone du système solaire, où notre soleil aurait le diamètre d’un ballon de Basket-Ball, la Terre serait à 27 mètres et aurait la taille d’une toute petite bille, de 2.3 millimètres de diamètre), on utilise des distances (et des volumes, qui proposer de représenter des distances par d’autres distances, bien plus petites.

Ici, l’exercice est, je pense, différent: on utilise le temps pour avoir une représentation d’une grandeur monétaire. Certes, un vieil adage prétend que le temps c’est de l’argent, mais je pense que la symbolique est vraiment forte ici. Si les décodeurs du Monde revenaient sur le calcul de cette grandeur (la “date de libération fiscale”) en le questionnant, personnellement c’est la symbolique attachée qui me fait m’interroger.

Les statistiques cachées derrière cette date du 29 juillet comme “date de libération fiscale”, ce sont

  • la dépense publique (ratio entre dépense publique et richesse créée) qui est de l’ordre de 57% en France (mais on ne voit pas vraiement le lien avec la fiscalité ici)
  • le “taux de taxation réel du salarié moyen” utilisé par certains, qui est aussi de l’ordre de 57% en France

Et 57% d’une année, c’est (en gros) entre le 1er janvier et fin juillet, d’où la date mentionnée dans l’article. Symboliquement, que lit-on, avec cette seconde statistique ? Qu’entre le 1er janvier et le 29 juillet on travaille “pour l’état”, pour payer ses impôts, et qu’ensuite, on peut (enfin) travailler “pour soi”. Au début, on travaille pour “les autres” par opposition à la fin où on travaille pour “soi”. Avec cette représentation calendaire, on se dit qu’on perd son temps, probablement. Mais à y regarder de plus près, est-ce que je travaille vraiment pour moi ensuite ? Personnellement, je pense que je passe aussi quelques jours à travailler pour ma banque, pour payer les intérêts du crédit que je rembourse tous les mois pour avoir acheté mon logement il y a quelques temps. Parce que mine de rien, des charges, il y en a beaucoup, il n’y a pas que l’Etat. Pour avoir cotisé pour ma retraite lorsque je travaillais au Canada, puis pour une assurance maladie, j’ai bien vu qu’à l’époque, je travaillais peut-être moins pour payer ces “charges fiscales”, mais je devais travailler plusieurs jours par an (voire semaines) pour payer l’assurance maladie qui couvrait la famille. Bref, le “pour moi” me dérange. Et que dire du “pour les autres” ? La base de la fiscalité, c’est la redistribution, non ? C’est la solidarité… Je ne paye pas des impôts pour les autres, mais “pour nous”, non ?

Bref, cette séparation binaire du temps me paraît fallacieuse et dangereuse car au lieu d’éclairer, elle n’a qu’une visée idéologique…. et c’est dommage.

Tests non-paramétriques et simulations

Lors du dernier cours de statistique, nous avons présenter les tests d’ajustment de lois. Nous avions illustré ces tests à partir de la taille d’individus (déjà utilisé pour présenter l’ajustement de lois, et l’estimation de densité), correspondant à https://latex.codecogs.com/gif.latex?\boldsymbol{x}=\{x_1,\cdots,x_n\}.

> Davis=read.table(
+ "http://socserv.socsci.mcmaster.ca/jfox/Books/Applied-Regression-2E/datasets/Davis.txt")
> Davis[12,c(2,3)]=Davis[12,c(3,2)]
> X=Davis$height
> n=length(X)

On notera https://latex.codecogs.com/gif.latex?%20(x_{i:n}) la statistique d’ordre, au sens où

https://latex.codecogs.com/gif.latex?%20x_{1:n}\leq%20x_{2:n}\leq\cdots\leq%20x_{n-1:n}\leq%20x_{n:n}

Parmi les outils graphiques, nous avons vu le PP plot (graphique probabilité-probabilité) et le QQ plot (graphique quantile). Le code pour créer un PP plot peut être le suivant

> PP=function(Y,F=pnorm){
+   n=length(Y)
+   x=F(sort(Y))
+   y=seq(1/n/2,1-1/n/2,by=1/n)
+   return(data.frame(x=x,y=y))
+ }

qui représente (à un détail près) le nuage de points

https://latex.codecogs.com/gif.latex?\left\{F_0(x_{i:n}),\frac{i}{n}\right\}

et celui pour le QQ plot serait

> QQ=function(Y,Q=qnorm){
+   n=length(Y)
+   x=Q(seq(1/n/2,1-1/n/2,by=1/n))
+   y=sort(Y)
+   return(data.frame(x=x,y=y))
+ }

qui représente (toujours à un détail près) le nuage de points

https://latex.codecogs.com/gif.latex?%20\left\{F_0^{-1}\left(\frac{i}{n}\right),x_{i:n}\right\}

où https://latex.codecogs.com/gif.latex?%20F_0 est la loi que l’on cherche à tester, au sens où on a https://latex.codecogs.com/gif.latex?%20H_0:F=F_0 avec comme hypothèse alternative https://latex.codecogs.com/gif.latex?%20H_1:F\neq%20F_0.

Continue reading Tests non-paramétriques et simulations

Autres Données pour le cours de Statistiques

Encore quelques données pour faire un TP de statistique. Ce sont des données de coûts d’ouragans, aux Etats-Unis,

> library(gdata)
> db=read.xls(
+     "http://sciencepolicy.colorado.edu/publications/special/public_data_may_2007.xls",
+     sheet=1)
> stupidcomma = function(x){
+  x=as.character(x)
+  for(i in 1:10){x=sub(",","",as.character(x))}
+  return(as.numeric(x))}
> base=db[,1:4]
> base$Base.Economic.Damage=
Vectorize(stupidcomma)(db$Base.Economic.Damage)
> base$Normalized.PL05=Vectorize(stupidcomma)(db$Normalized.PL05)
> base$Normalized.CL05=Vectorize(stupidcomma)(db$Normalized.CL05)

Lire la base est un peu pénible, à cause de virgules qui traînent comme séparateur de milliers. Mais peu importe, le petit code permet de nettoyer la base.

  • Fréquence d’ouragans

La première étape pourra être de travailler sur la fréquence. On pourra ajuster des lois de Poisson, ou binomiales, ou binomiales négatives, et faire des tests de stabilité du paramètre, avec le temps (on pourra  couper en 4 ou 5 blocs).

> TB <- table(base$Year)
> years <- as.numeric(names(TB))
> counts <- as.numeric(TB)
> years0=(1900:2005)[which(!(1900:2005)%in%years)]
> db <- data.frame(years=c(years,years0),
+       counts=c(counts,rep(0,length(years0))))
> db[88:93,]
   years counts
88  2003      3
89  2004      6
90  2005      6
91  1902      0
92  1905      0
93  1907      0
> plot(years,counts,type='h')

  • Coût des ouragans

Dans un second temps, on pourra travailler sur le coût des ouragans

> plot(base$Normalized.PL05/1e9,type="h")

On pourra tenter d’ajuster des lois, log normale, gamma, Pareto, etc. Faire des tests d’ajustement. Et tester la stabilité des paramètres, dans le temps (en coupant en 2 blocs, par exemple, ou davantage).

Données pour le cours de Statistique

Pour les derniers cours de Statistique, j’avais envoyé  une feuille avec quelques exercices, à regarder. J’en profile pour mettre en ligne les exercices sur l’inférence, et ceux sur les tests. Je mets aussi en ligne un tableau contenant les résultats de cinq établissements sur Rennes, au Brevet des Collèges. On pourra utiliser ces données pour faire un peu de statistique.

Test du chi-deux et indépendance

Considérons le tableau de contingence suivant

> N=margin.table(HairEyeColor, c(1,2))
> N
       Eye
Hair    Brown Blue Hazel Green
  Black    68   20    15     5
  Brown   119   84    54    29
  Red      26   17    14    14
  Blond     7   94    10    16

avec ici les comptages, que l’on peut aussi traduire sous la forme de probabilités jointes (empiriques)

> N/n
       Eye
Hair    Brown  Blue Hazel Green
  Black 0.115 0.034 0.025 0.008
  Brown 0.201 0.142 0.091 0.049
  Red   0.044 0.029 0.024 0.024
  Blond 0.012 0.159 0.017 0.027

Continue reading Test du chi-deux et indépendance

Rappels de statistique

Jeudi et vendredi, on finira les rappels de cours, avec quelques rappels de statistique. Les slides sont en ligne,

Sinon, histoire d’illustrer un peu les différents aspects, j’avais dit que j’aimais bien utiliser des simulations. En particulier, il sera intéressant de générer plusieurs (beaucoup) échantillons de taille donnée. Le code R est le suivant (une ligne correspond à un échantillon, de 20 observations)

> n=20
> m=1000
> M=matrix(rnorm(n*m),m,n)

On utilise ensuite la statistique que l’on souhaite étudier, qui peut être la moyenne empirique,

> Xbarn=apply(M,1,mean)

ou la variance

> Sn=apply(M,1,var)

On peut ensuite visualiser la distribution de la statistique, sur tous les échantillons. Par exemple, le fait que l’estimateur de la moyenne est gaussien, on peut faire un histogramme

> hist(Xbarn,probability=TRUE,
+ breaks=seq(-1,1,by=.1),col="light blue")

et rajouter la densité de la loi normale (car la théorie nous dit que la moyenne suit une loi normale),

> u=seq(-2,2,by=.01)
> lines(u,dnorm(u,sd=1/sqrt(n)),col="red",lwd=2)

Cela dit, comme Ewen Gallic fera essentiellement du ggplot2 dans son cours de R, je vais aussi mettre les codes en ggplot2 (et essayer de m’y convertir… ce n’est pas l’envie qui manque !)

> ggplot(data=data.frame(Xbarn = Xbarn), aes(x = Xbarn)) +
+   geom_histogram(aes(y = ..density..),
+   binwidth = .1, fill = "light blue", col = "black") +
+   stat_function(fun=function(x) dnorm(x,sd=1/sqrt(n)),
+ col = "red", size = 1) + ylab("Densite") +
+   ggtitle("Histogramme de Xbarn")

Cela dit,  au lieu de représenter la loi normale, on peut aussi visualiser un estimateur non paramétrique de la densité (on en reparlera dans quelques semaines), comme l’estimateur à noyau

> ggplot(data=data.frame(Xbarn = Xbarn), aes(x = Xbarn)) +
+   geom_histogram(aes(y = ..density..),
+   binwidth = .1, fill = "light blue", col = "black") +
+   geom_line(stat="density", col = "red", size = 1) +
+   ylab("Densite") +
+   ggtitle("Histogramme de Xbarn")

On peut aussi voir que la variance empirique suit une loi du chi-deux,

> hist((n-1)*Sn,probability=TRUE,
+ breaks=seq(0,56,by=2),
+ col="light blue")
> u=seq(0,56,by=.1)
> lines(u,dchisq(u,df=n-1),col="red",lwd=2)

avec là encore une version en ggplot2

> ggplot(data = data.frame(x = (n-1)*Sn), aes(x = x)) +
+   geom_histogram(aes(y = ..density..),
+   binwidth = 2, fill = "light blue", col = "black") +
+   geom_line(stat="density", col = "red", size = 1) +
+   xlab(expression((n-1) %*% Sn)) +
+   ylab("Densite") +
+   ggtitle(expression(paste("Histogramme de ", (n-1) %*% Sn)))

On peut aussi utiliser ces simulations pour comprendre l’idée d’intervalle de confiance à 90%,

> x1=Xbarn-1.64*1/sqrt(n)
> x2=Xbarn+1.64*1/sqrt(n)
> ind=x1*x2>0
> plot(1:m,Xbarn,ylim=c(-1,1),col="white",
+ ylab="intervalle de confiance",xlab="")
> segments(1:m,x1,1:m,x2,col=c("light blue","red")[ind+1],lwd=1.5)
> abline(h=0,lty=2)

ou pour la version ggplot2 (mais on reparlera de tout ça en cours, pas forcément cette semaine cela dit)

> ggplot(data = df) +
+   # Les segments. La couleur doit etre fonction d'un logique
+   geom_segment(aes(x = x, xend = x, y = x1, yend = x2, col = x1*x2>0)) +
+   # La ligne horizontale en y = 0
+   geom_hline(yintercept = 0, linetype = "dashed") +
+   # Changer les couleurs des segments
+   scale_colour_manual(values = c("light blue" ,"red")) +
+   # Modifier le titre de l'axe des y
+   ylab("Intervalle de confiance") +
+   # Modifier le titre de l'axe des x
+   xlab("") +
+   # Cacher la legende
+   theme(legend.position="none")

On note que dans 90% des scénarios (notion d’intervalle de confiance à 90%), la vraie valeur (connue ici car on simule nos échantillons) est effectivement dans l’intervalle de confiance.

Faire parler les chiffres… n’importe comment

Cette fin de semaine, Martin Grandjean a mis en ligne un billet intéressant sur son blog, sur l’utilisation des statistiques (dans un but de propagande). L’exercice n’est pas nouveau, mais Martin soulève des questions, malheureusement importantes et complexes. Dans un paragraphe, intitulé “faire parler les chiffres… n’importe comment” (que j’ai repris comme titre, j’avoue avoir hésité avec “with great power comes great responsibility“), on retrouve l’analyse (rapide) d’un graphique, présenté ci-dessous. Le graphique parle de “statistiques” liées à des problèmes démographiques et d’immigration, sauf que comme le note Martin, il y a un “petit problème” : “ces statistiques n’en sont pas vraiment. Ou plutôt, il s’agit d’une extrapolation libre et “linéaire”“.

Je reviendrais sur ces points dans deux minutes, avant je voulais revenir sur une petite phrase que l’on retrouve un peu plus loin: “comment peut-on se permettre de prédire une tendance sur cinquante ans à partir d’une tendance observée sur dix ans ?“. Cette remarque est importante… car beaucoup de monde se pose cette question, et les statisticiens se doivent d’apporter une réponse. Ce problème m’avait été posé (plus ou moins) par une journaliste il y a quelques années maintenant, sur la solvabilité des compagnies d’assurance. Avec Solvabilité II, en Europe, les compagnies doivent calculer des VaR à 99.5%, c’est à dire liés à des événements ayant une période de retour de 200 ans. Avec au mieux 25 ans de données ! Dans ce cas, la nuance vient du fait que l’on cherche à estimer un quantile associé à une probabilité faible avec peu de données (et que l’on assimile rareté et période temporelle). C’est le même problème qui se pose en hydrologie, quand on veut construire une digue suffisamment haute pour que seul un événement millénaire provoque une inondation, mais que l’on doit prévoir avec 50 ans de données. Cela dit, quand un assureur fait de la mortalité prospective, il a souvent 20 ou 30 ans de données, et doit estimer la probabilité qu’un assuré de 25 ans soit encore en vie dans 50 voire 75 ans. C’est le genre de message que j’essayais de faire passer il y a quelques années dans une table ronde sur la pérennité des régimes de retraites, en France. Bref, cet exercice est un vrai exercice. Et il convient de bien comprendre le modèle pour comprendre ce qu’on fait ! Dire que c’est stupide ne pourra pas suffire.

Car quand on a peu de données (j’ai envie de dire que c’est le cas ici), c’est le modèle qui va imposer les conclusions: tous les modèles semblent valides quand on a peu de données, et les conclusions dépendront plus du modèle retenu que des données. C’est moins le cas quand on a des gros volumes de données (sans parler de big-data, je garde ça pour un autre billet).

Revenons sur le problème en question, car il est intéressant. On nous parle ici d'”extrapolation linéaire” dans la publicité. “Extrapolation” signifie que l’on va construire un modèle, et qu’on va l’utiliser pour faire une prévision. Et il y a “linéaire“. A priori, un modèle linéaire, ce n’est pas trop ambigu… sauf qu’on peut faire un modèle sur le niveau, ou sur le taux de croissance. Si on fait un modèle très simple – constant – sur le taux de croissance, on aura alors une croissance exponentielle, c’est aussi simple que ça. Donc le modèle peut être linéaire, ou exponentiel (on regardera un peu les deux).

Autre point. On a ici trois populations,

  • les suisses, 
  • les étrangers, 
  • la population totale, 

Le graphique est un peu trompeur… Il donne l’illusion que l’on extrapole la première variable  et la troisième, . Le danger est qu’il convient, dans ce cas, de contraindre un minimum les modèles, car . On pourrait très bien imaginer qu’en extrapolant, on ait . Donc il faudra faire un peu attention… Le plus simple sera probablement de construire des modèles sur  et  en s’assurant de la positivité de ces deux variables.

Dernier point. A-t-on le droit de faire des modèles indépendants. Peut-on supposer que les deux séries  et  évoluent indépendamment l’une de l’autre ? A priori non, et dans ce cas, c’est un peu plus complexe, car il faut un modèle bivarié. Bref, je trouve que ce petit exemple donne matière à réflexion sur ce qu’est la modélisation statistique. Et sur l’incroyable pouvoir qu’on les statisticiens et les économètres (car ceux sont eux qui construisent le modèle… “with great power comes great responsibility” comme dirait l’autre).

Car il faut que les choses soient claires. Il est impossible d’être neutre ! Ça n’a pas de sens… On peut construire plusieurs modèles, qui sembleront tous valides, et ensuite utiliser un critère de choix de modèles (les plus usuels étant peut-être le critère AIC d’Akaike, ou le BIC de Schwarz). Mais là aussi, je n’ai aucun moyen pour dire : ce modèle est le meilleur, et il est neutre. Non, il sera le meilleur pour un critère de choix, que je choisis.

Histoire d’être plus clair, essayons de faire des régressions. Le modèle le plus simple serait un modèle linéaire,

avec des résidus – a priori – non corrélés, . C’est que donnerait deux régressions linéaires sous Excel, pour ceux qui ne sont pas à l’aise avec la modélisation. Classiquement, ce modèle s’écrit

avec  (conditionnellement à notre variable explicative). Lançons nous… Bon, quand j’ai cherché les données, je n’ai pas trouvé 3 observations, mais une bonne trentaine pour la population résidente, et les étrangers. On va donc utiliser toutes les données pour faire notre projection (et éventuellement voir ce que donnerait une régression avec seulement 3 observations)

> D=read.table("http://freakonometrics.free.fr/suisse.csv",sep=",",
+ header=TRUE)
> tail(D)
      X      N1      N3      N2
27 2006 1554527 7508739 5954212
28 2007 1602093 7593494 5991401
29 2008 1669715 7701856 6032141
30 2009 1714004 7785806 6071802
31 2010 1766277 7870134 6103857
32 2011 1815994 7954662 6138668

Comme on a un modèle linéaire, on peut faire indifféremment l’estimation du modèle sur deux des trois séries. Par simplicité, modélisons des deux séries qui nous intéressent. L’estimation donne ici

> reg2=lm((N2/100000)~X,data=D)
> reg3=lm((N3/100000)~X,data=D)
> pred2=predict(reg2,newdata=data.frame(X=1980:2060))/100000
> pred3=predict(reg3,newdata=data.frame(X=1980:2060))/100000
> lines(1980:2060,pred2,col=COL[1],lwd=2)
> lines(1980:2060,pred3,col=COL[6],lwd=2)

Maintenant, il faut s’interroger sur la pertinence de notre modèle. Pour ma part, j’ai rarement vu des modèles linéaires pour modéliser des des évolutions de populations. Mais sur du court terme (une centaine d’années), pour quoi pas. Même avec un modèle de type Maltusien, c’est possible. Pour ma part, quand on travaille sur des données de comptage, j’aime beaucoup la régression de Poisson,

A la différence de notre premier modèle, on aura ici de l’hétéroscédasticité : plus la population grande, plus la variance (et donc le terme d’erreur) sera grande. C’est plus réaliste que le premier modèle. Rien de plus facile que d’estimer ce modèle

> reg2=glm((N2/100000)~X,data=D,family=poisson(link="identity"))
> reg1=glm((N1/100000)~X,data=D,family=poisson(link="identity"))

On peut en profiter d’ailleurs pour regarder un peu les intervalles de confiance (sur notre modèle) et relativiser un peu les conclusions

> pred2p=predict(reg2,newdata=data.frame(X=1980:2060),type="response",se.fit=TRUE)
> pred1p=predict(reg1,newdata=data.frame(X=1980:2060),type="response",se.fit=TRUE)

La taille de l’intervalle de confiance nous suggère ici d’être prudent avant d’affirmer quoi que ce soit… et encore, on a de la chance, on a ici 35 ans de données et pas 3 !

> I=which(D$X%in%c(1990,2000,2010))
> reg2=glm((N2/100000)~X,data=D[I,],family=poisson(link="identity"))
> reg1=glm((N1/100000)~X,data=D[I,],family=poisson(link="identity"))

Bon, classiquement, quand on fait une régression de Poisson, on utilise un lien exponentiel,

On va alors s’assurer que les deux comptages soient toujours positifs ! Ce qui est bien compte tenu de notre problème. Mais en contrepartie, on a la contrainte d’avoir une croissance (ou une décroissance) exponentielle de nos populations. On n’a pas trop le choix (car quand on y pense, avec un modèle linéaire, il y a forcément une date, avant laquelle, ou après laquelle, la population sera négative… ce qui est gênant)

> reg2=glm((N2/100000)~X,data=D,family=poisson)
> reg1=glm((N1/100000)~X,data=D,family=poisson)
> pred2p=predict(reg2,newdata=data.frame(X=1980:2060),type="response")
> pred1p=predict(reg1,newdata=data.frame(X=1980:2060),type="response")

(la courbe en pointillé en arrière plan est le modèle linéaire).

On continue ? On peut en fait avoir un modèle bivarié, pour dire que nos deux courbes évoluent ensemble. Il existe par exemple le modèle de Poisson à choc commun, de densité jointe

avec . Sous R, le package qui permettait de faire cette régression n’est plus disponible, mais peu importe, on peut le recoder, ça prendra trois minutes,

> f=function(Z,L){
+ si=function(i) choose(Z[1],i)*choose(Z[2],i)*gamma(i+1)*
+ (L[3]/(L[1]*L[2]))^i
+ s=Vectorize(si)(0:min(Z))
+ p=exp(-sum(L))*L[1]^Z[1]/gamma(Z[1]+1)*L[2]^Z[2]/gamma(Z[2]+1)*sum(s)
+ return(p)
+ }
> minuslogL=function(B){
+ h=function(x) exp(x)
+ logL=function(i) log(f(Y[i,],
+ c(h(B[1]+B[2]*X[i]),h(B[3]+B[4]*X[i]),h(B[5]+B[6]*X[i]))))
+ return(-sum(Vectorize(logL)(1:n)))
+ }
> optim(c(lm(log(Y[,1])~X)$coefficients,
+ lm(log(Y[,2])~X)$coefficients,0,0),minuslogL)->maxL
> Bstar=maxL$par
> Bstar
  (Intercept)             X   (Intercept)             X               
-4.343823e+01  2.303903e-02 -3.430377e+00  3.746384e-03 

 1.506016e-02 -9.743153e-04 
> predbiv=function(x,B=Bstar){
+ h=function(x) exp(x)
+ return(c(h(B[1]+B[2]*x)+h(B[5]+B[6]*x),
+          h(B[3]+B[4]*x)+h(B[5]+B[6]*x)))

Si on construit notre prévision, on est ici proche de notre précédant modèle,

On pourrait aussi utiliser un modèle de série temporelle: dire que l’on a une tendance linéaire, ou exponentielle, et dire que le bruit peut être modélisé par un modèle autorégressif (bivarié ou pas).

avec un modèle autorégressive pour le bruit,

ou de manière beaucoup plus générale,


i.e. un modèle autorégressif vectoriel,

On peut estimer ce modèle assez facilement,

> reg2=glm(N2~X,data=D,family=poisson)
> reg1=glm(N1~X,data=D,family=poisson)
> Z=cbind(residuals(reg2),residuals(reg1))
> library(vars)
> regvar=VAR(Z,p=1)

Si on regarde le modèle, on notera que la matrice d’autocorrélation est bien pleine, avec de la causalité dans tous les sens (ce qui devrait augmenter encore la taille de nos intervalles de confiance)

> summary(regvar)

VAR Estimation Results:

Estimation results for equation y1: 
=================================== 
y1 = y1.l1 + y2.l1 + const 

      Estimate Std. Error t value Pr(>|t|)    
y1.l1  1.13635    0.05191  21.891   <2e-16 ***
y2.l1 -0.11922    0.03411  -3.495   0.0016 ** 
const  0.73037    0.52316   1.396   0.1737    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Estimation results for equation y2: 
=================================== 
y2 = y1.l1 + y2.l1 + const 

      Estimate Std. Error t value Pr(>|t|)    
y1.l1  0.44698    0.12384   3.609  0.00118 ** 
y2.l1  0.83294    0.08139  10.234 5.76e-11 ***
const  0.87575    1.24811   0.702  0.48868    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Je continue ? Allez, un petit dernier exemple pour la route ! Regardons le petit modèle suivant,

Avec ce modèle, le ratio entre les deux populations sera toujours le même. Ce qui fera taire tous les débats… et pour l’estimer, il suffit de rajouter la contrainte dans les codes précédents (avec l’estimation par maximum de vraisemblance)

Maintenant, comme on avait des séries temporelles, je n’ai pas présenté de modèles nonparamétriques, car c’est délicat de faire de la prévision avec ces modèles. Par contre, si on avait des données individuelles, on aura pu s’en donner à cœur joie, à faire du lissage dans tous les sens ! Ce qui nous aurait donné encore plus de latitude dans notre modèle. Je passe aussi sous silence ici que mes modèles sont probablement trop simples. Je veux dire par là – je connais mal les procédures pour avoir la nationalité suisse – que quand ses enfants, et ses petits-enfants, sont nés en Suisse, j’ai du mal à me dire que la personne est un “étranger“. Au bout de quelques années, on peut imaginer un modèle dynamique où les étrangers passent de “étranger” à “Suisse“. Mais j’ai peut-être un biais canadien dans ma perception de la nationalité.

Ce que je voulais dire, c’est que faire des prévisions à partir d’un (relativement) petit jeu de données n’est jamais simple, et que le choix du modèle va complètement déterminer la visualisation que l’on veut avoir pour nos données, et du type de prévision que l’on cherche. Et ne me demandez pas d’être neutre, c’est impossible… Une fois qu’on a compris ce problème, on comprend que certains chercheurs avancent que plus de la moitié des études publiées (souvent sur des petits échantillons) sont fausses. Je ne dirais pas pour ma part qu’elles sont “fausses“. Juste qu’il est facile, quand on a peu de données, de faire passer le message que l’on veut ! Ici, c’est facile à démonter (Martin l’a fait en quelques lignes, et il m’en aura fallu un peu plus car j’ai besoin d’un peu plus de temps – et de formalisme – pour raconter mes histoires), et il s’agissait uniquement d’un dessin de propagande. Il faut garder en mémoire que beaucoup de décisions de santé publique sont basés sur ce genre d’exercice. Avec en plus des données rarement publiques, et entachées de beaucoup de bruit… On pensera aux études sur les OGM, sur les antennes relais, sur le tabagisme (évoqué dans un précédant billet), etc. Et pour démonter ce genre d’étude, cela prend beaucoup plus de temps !

Les abaques comme outils de pédagogie statistique

Comme je l’avais mentionné à propos d’un vieux papier paru dans Le Canard, c’est souvent compliqué de communiquer sur les retraites. C’est un sujet complexe et assez rapidement technique. Historiquement, les abaques sont les tables de calculs,

Mais il faut admettre que retrouve un peu n’importe quoi dans ce terme d’abaque…. Certaines sont évoquées ici ou , voire pour quelque chose de très joli. On notera qu’il existe une idée sous-jacente commune, à savoir construire un outil graphique permettant de visualiser comment 3 paramètres peuvent être liés (généralement, ils sont liés par une fonction complexe).

  • Abaques et retraites

Parmi les outils que l’on retrouve dans les documents du COR, j’ai été surpris de voir ressortir l’utilisation des abaques (par exemple ici). On nous en fait même la promotion en vidéo,

Pour la petite histoire, j’avais évoqué son utilisation dans un graphique ici sur les problèmes d’échantillonnage. C’est effectivement un outil graphique intéressant, à condition de savoir le décoder…
L’idée est simple, à savoir lier les trois paramètres essentiels(ou supposés comme tels) du bilan d’un système de retraite, à savoir

  • le niveau des prélèvements,
  • le niveau de l’âge effectif de départ en retraite
  • le niveau des pensions.

La difficulté est de trouver une forme graphique permettant de relier ces trois niveaux (la dimension 3 étant souvent trop abstraite). Par exemple dans les sondages, pour construire un intervalle de confiance, les abaques sont utilisées (comme ici) pour relier

  • le niveau des intervalles de confiance (5%-95% ou 10%-90%)
  • le niveau de probabilité attendu dans le sondage
  • la taille de l’échantillon interrogé

On peut d’ailleurs trouver d’autres exemple en statistiques ici ou . Mais pour revenir au problème des retraite, sur le graphique ci-dessous (appendix 1 ici), en abscisse on met le taux de remplacement (c’est à dire le ratio de la pension moyenne par rapport au revenu moyen) et en ordonnée on met le niveau d’un hausse du taux de cotisation (en points de cotisation). chaque ligne représente un âge moyen de départ en retraite effectif.

On se fixe alors des scénarios d’équilibre et on regarde la manière dont se déplacent les équilibres. J’avais prévu toute une analyse technique, mais je viens de me rendre compte qu’Antoine a mis en ligne un billet remarquable sur le sujet sur http://www.ecopublix.eu/2010/…, que je ne pourrais jamais égaler en clarté. Donc je vais proposer une autre construction d’abaque afin de montrer son intérêt pédagogique.

  • Comprendre l’actualisation

Sur le site http://images.math.cnrs.fr/…, Xavier Caruso explique son émerveillement devant la difficulté de faire des calculs d’actualisation. En un sens, ça me rassure de voir que les calculs actuariels ne paraissent pas triviaux à tout le monde (et surtout plein de contre-intuitions). Xavier essaye d’expliquer la différence entre

  • ouvrir un compte rémunéré à 4,5% sur lequel on dépose 500 euros chaque mois pendant 10 ans
  • ouvrir un compte rémunéré à 4,5% sur lequel on dépose 60000 euros initialement sans plus jamais l’alimenter ensuite.

Si effectivement la somme des 120 versements de 500 euros correspond aux 60000 euros, la valeur intégrant les intérêts n’est pas du tout la même au bout de 10 ans, car dans le premier car, comme le disent les banquiers, “les intérêts se capitalisent“. Bon, si Xavier fait des maths (et donc des développements limités par avoir des formules fermés), on va plutôt faire ici un peu de calculs numériques. Tout d’abord, pour visualiser l’évolution des montants disponibles sur les comptes, rien de plus simple,

> tps=1/12*(1:120)
> plot(tps,60000*(1+0.045)^tps,ylim=c(0,100000))
> lines(tps,cumsum(500*(1+0.045)^(tps)))

les montants au bout de 10 ans étant

> sum(500*(1+0.08525)^(1/12*(1:120)))
[1] 93180.33
> sum(500*(1+0.045)^(1/12*(1:120)))
[1] 75514.32

Ce qui confirme l’intuition que nous avions. Comme le note Xavier, il faudrait un taux d’intérêt deux fois plus grand pour le placement mensuel, la preuve

> sum(500*(1+0.09)^(1/12*(1:120)))
[1] 95543
> sum(500*(1+0.08525)^(1/12*(1:120)))
[1] 93180.33

Xavier obtient cette valeur en effectuant des développements limités. Il résume ça à l’aide du dessin

On peut aller plus loin, en notant qu’à taux sur le placement mensuel donné, ainsi que la maturité, le taux sur le placement bloqué qui rapporte la même valeur à échéance est donné par

 

>1))/(500*12*10))^(1/10)-1
[1] 0.04762238

(dans le cas d’un placement rapportant 9%). En fait, si on regarde le montant dont on dispose sur les comptes à chaque date, on obtient précisément les valeurs suivantes, en fonction du temps,

on retrouve effectivement qu’à 4,5%, les deux placements ne donnent pas le même montant à échéance (ici 10 ans). Alors que si le placement alimenté tous les mois rapportait 8,525%, on aurait exactement la même somme au bout de 10 ans

ou si on compare avec un placement à 9% (le double, ce qui est obtenu par développement limité), on obtient

Pour trouver la valeur exacte du multiplicateur en fonction du taux offert sur le compte bloqué et la maturité, la fonction est tout simplement

> taux=4.5/100
> T=10
> f=function(k){(1*12*T)*(1+taux)^T-sum(1*(1+k*taux)^(1/12*(1:(12*T)))) }
> uniroot(f,interval=c(0,20))$root
[1] 1.894332

L’idée d’utiliser une abaque peut s’avérer intéressante, car on essaye ici de relier

  • le taux de rendement de l’argent bloqué
  • le taux proposé pour le placement alimenté tous les mois
  • la durée envisagée du placement

Je peux mettre en abscisse le taux du placement bloqué. En ordonné, je met le multiplicateur utilisé pour le placement alimenté tous les mois. A maturité donné (par exemple 10 ans en rouge), on peut regarder la valeur du multiplicateur pour que les valeurs à échéance coïncident,

On retrouve directement que le taux équivalent pour un placement alimenté tous les mois doit être 1,89 fois le taux offert sur le placement bloqué si on se fixe une échéance de 10 ans pour égaliser les placements. Le code est simplement le suivant

> tx=seq(0.005,.05,by=.001)
> tps=seq(5,25,by=.25)
> M=matrix(NA,length(tps),length(tx))
> for(i in 1:length(tps)){
+ for(j in 1:length(tx)){
+ T=tps[i]
+ taux=tx[j]
+ f=function(k){
+ (1*12*T)*(1+taux)^T-sum(1*(1+k*taux)^(1/12*(1:(12*T))))
+ }
+ M[i,j]=uniroot(f,interval=c(0,2))$root
+ }}
> contour(tps,tx,M,lwd=2,
+ xlab="Durée du placement (en années)",ylab="Taux (du placement bloqué)")
> abline(v=seq(5,25,by=1),lty=2,col="grey")
> abline(h=seq(0,.05,by=.005),lty=2,col="grey")

Bref, l’actualisation est un sujet délicat (même ce n’est que des calculs de sommes et de puissances). Et faire des petits dessins permet souvent de mieux comprendre…

  1. sum(500*(1+0.09)^(1/12*(1:120 []

Discussions géométriques du matin (partie 2)

Suite rapide de mon billet sur la construction de triangles avec des bouts de bois (ici). Adrian me faisait remarquer () qu’il y avait plus simple pour calculer la probabilité de construire un triangle (je suis d’ailleurs affreusement jaloux de ne pas y avoir pensé tout seul). Mon principal soucis que j’ai en lisant le commentaire est “ça semble juste, mais est-ce que c’est juste ?” (mon expérience en TD de proba est que l’on peut souvent trouver un paquet de solutions qui “semblent” justes). Bref, faisons preuve de rigueur quelques minutes (même si France Inter continue à nous dire que le mot est toujours tabou ces temps-ci).
Rappelons que si on dispose d’un échantillon L'image “https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-01.png” ne peut être affichée car elle contient des erreurs. i.i.d. de loi https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-02.png, et de densité https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-03.png, et que l’on définie les statistiques d’ordre

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-04.png

la loi du couple https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-05.png, avec https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-06.png,est alors

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-07.png

à condition que https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-08.png (sinon la densité sera nulle). C’est jsute des histoires de compter les classements possibles, et puis de calculer quelques probas simples. Un cas particulier connu de tous est la loi du couple https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-09.png, i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-10.png

Dans le cas qui nous intéresse, supposons que l’on tire deux bâtons, alors la loi du couple s’écrit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-11.png

alors que la densité du maximum était

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-12.png

donc la loi du minimum conditionnellement au maximum doit d’écrire, par la formule de Bayes

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-13.png

qui est effectivement la densité de la loi uniforme sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-14.png. Bref, ça marche pour deux bâtons.
Dans le cas de trois bâtons, la loi du triplet s’écrit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-15.png

alors que la loi du maximum devient

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-30.png

aussi la loi conditionnelle du couple des deux premiers s’écrit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-31.png

Or si l’on avait deux lois uniformes sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-16.png, compte tenu de ce que l’on a écrit auparavant, la loi du couple des deux statistiques d’ordre seraient

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-17.png

soit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-18.png

Bref, il semble qu’effectivement, le résultat soit correct: conditionnellement au maximum, les autres variables sont indépendantes, et suivent un loi uniforme sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-19.png.
Moralité, effectivement, les non-maximum sont indépendants, et leur loi est uniforme sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-16.png. On peut alors sommer ces deux lois uniformes et indépendantes, ce qui donne une loi triangulaire, sur https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-20.png, avec le sommet en https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-21.png. Et par des raisons de symmétrie, la probabilité de dépasser https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-21.png est alors https://perso.univ-rennes1.fr/arthur.charpentier/latex/tr-ordre-22.png.Merci Adrian. C’est effectivement plus élégant que mon calcul moche…

Estimation de quantile par noyau beta

Le papier sur l’estimation de quantile par noyau beta, coécrit avec Abder Oulidi, est accepté pour publication dans Statistics and Computing, http://link.springer.com/…

In this paper we propose several nonparametric estimators of quantiles based on Beta kernel and applied to transformed data by the generalized Champernowne distribution initially fitted to the data. A Monte-Carlo based study will show that those estimators improve the efficiency of a traditional ones, not only for light tailed distributions, but also heavy tails, when the probability level is close to 1. We also compare these estimators with the Extreme Value Theory Quantile applying to Danish data on large fire insurance losses.