Category Archives: Statistics

Prédire ou Prévoir, réflexions (probablement) personnelles

Il y a quelques années, j’avais fait un billet où je m’énervais contre une phrase que j’avais entendu (“l’homme qui avait prédit la crise). Le point de départ était surtout un énervement lié au fait que beaucoup de journalistes vivent avec des oeillères (par exemple en associant “économiste” à “type à la télé qui parle d’économie“), et ne connaissent pas la diversité de la communauté des économistes (surtout que la pluralité n’est pas vraiment une vertue présente dans les grands médias, la manière dont sont présentés les différents candidats à la future élection présidentielle peut en être une illustration). Mais je m’écarte…. en fait, un de mes points portait sur l’utilisation du terme “prédit“.

Avec un peu de recul, je ne sais pas si la langue française à changé depuis, mais j’aurais envie de revenir sur ma distinction entre “prévoir” et “prédire“. Prévoir et prédire ont deux racines latines distinctes, prévoir venant du latin praevidere, et prédire de praedicere. Beaucoup notent que “prédire” a longtemps eu le sens de “prophétiser“, où l’annonce (de ce qui doit arriver) est souvent d’inspiration divine. En anglais, c’est presque plus compliqué (car il y a plus de synonymes) mais la même confusion existe. On peut retrouver “foretell” construit sur “fore” et “tell” (qui est à relier à une notion de prémonition), “forespell” construit sur “fore” et “spell” (on retroue la notion de sort, proche du néerlandais “voorspellen”), “foresee” construit sur “fore” et “see”, ou encore  “forecast“, dont on reparlera très bientôt. A côté de ces mots, construits sur le préfixe “fore”, il existe le mot “predict“, tiré du latin praedicere, construit sur “prae” ‎(avant, que l’on retrouvé dans “pré” en français) et “dicere” ‎(dire). Comme le note le dictionnaire, “predicted” a un sens proche de “foretold” mais “without the mystical flavour”. Bref, “prédit” a longtemps eu une connotation mystique, divine.

Pourtant, les deux semblent parfois être interchangeables. Au point que fin 2006, Le Devoir écrivait « parce que cette administration n’a pas su prévoir ce qui était prévisible, nous en sommes rendus à siéger à 6 h un lundi matin! ». “prévoir ce qui était prévisible” vaste programme, non ?

Si on reste sur l’anglais, on peut tenter une distinction. Par exemple “forecast“, c’est à dire “prévoir“, est associé aux ouvrages suivants, sur amazon.com

Autrement dit, “forecast” (ou “prévision“) correspond à un aspect temporel. On cherche à prévoir, à un instant donné https://latex.codecogs.com/gif.latex?t la valeur d’une variable https://latex.codecogs.com/gif.latex?x à une date ultérieur, disons https://latex.codecogs.com/gif.latex?t+h. Classiquement, on cherche (et ça sera notre “prévision”)

https://latex.codecogs.com/gif.latex?\mathbb{E}_t[X_{t+h}]=\mathbb{E}[X_{t+h}|x_t,x_{t-1},\cdots]

C’est pour cela que bon nombre de cours de séries temporelles, ou séries chronologiques sont des cours de “modèles de prévision“.

Mais on commence à voir arriver des ouvrages de “méthodes prédictives“, “predictive modeling“, proposant des “prédictions

Si on regarde de plus près, ces techniques sont utilisées non pas dans un contexte de séries temporelles, mais de données individuelles. Ici, on cherche à prévoir une variable d’intérêt https://latex.codecogs.com/gif.latex?y à l’aide de covariables (variables dites explicatives https://latex.codecogs.com/gif.latex?x_1,\cdots,x_k. La prédiction sera classiquement

https://latex.codecogs.com/gif.latex?\mathbb{E}_{\boldsymbol{x}}[Y]=\mathbb{E}[Y|x_1,\cdots,x_k]

Sur le graphique du taux d’imposition, en France

Voilà plusieurs années (en 2011) Landais, Piketty et Saez publiaient le graphique suivant, présentant le sytème fiscal, en France, avec l’évolution du taux global d’imposition, en fonction du revenu.

https://freakonometrics.hypotheses.org/files/2017/02/revol-fiscale.png

Il y a deux ans, Yves Guiard critiquait ce graphique (rappelant que “la conversion numérique/graphique obéit à certaines règles“) et proposait une version “correcte“, sous la forme

J’ai déjà dit – voilà plusieurs semaines – que cette vision dogmatique des règles de représentation graphique ne me plaisait pas. Mais le plus simple est peut être de reprendre la réponse de Thomas Piketty, sur cette critique (je me retrouve totalement dans le propos),

(je suis désolé de découvrir cette discussion aussi tard, Béatrice m’a posté un commentaire ce matin, qui pointait vers ces deux articles)

Mais essayons de regarder un peu plus en détails. Oui, si on utilise la représentation standard de la courbe, on retrouve le graphique suivant

avec les plus pauvres à gauche, et les plus riches à droites. On a ici un graphique où le même pas sur l’axe des abscisses correspond au même nombre d’individus. Les 10% les plus pauvres (bande verticale la plus à gauche) ont un taux moyen d’imposition de 40.46% alors qu’il est de 48.09% pour les 10% qui sont juste au dessus du salaire médian (bande juste à droite de 50%). Je ne reviendrais pas sur le débat de commencer le graphique à 0% (ce débat aurait du sens, j’en conviens, mais commencer à 0% et mettre du texte entre l’axe des abscisses et la courbe, c’est un peu tricher), mais si on respecte cette règle, on obtient

Mais ce n’est pas cette visualisation qui est retenue dans le livre de Landais, Piketty et Saez, mais plutôt quelque chose qui ressemble au graphique ci-dessous (je commence à 30% pour éviter d’avoir un billet trop long)

Ce graphique est (je pense) plus intéressant que le premier. Bon, bien sûr, il se rapproche de la courbe parfaite que tout statisticien qui se respecte essaye reproduire

mais surtout, il permet de tenir compte de l’hétérogénéité qui existe dans ce dernier décile. Histoire d’illustrer mon propos, je voudrais proposer une représentation alternative de celle de Landais, Piketty et Saez. Pour se faire, on va utiliser une information complémentaire (en ligne sur le site) à savoir le revenu moyen, par centile, dans la population (voire plus fin dans le centile le plus élevé)

L’idée sera d’utiliser en abscisse non pas un pourcentage de personnes, mais un pourcentage de richesse (un peu comme avec la courbe de Lorenz),

Par exemple sur le graphique ci dessus, la bande la plus à gauche ne représente pas 10% des français, mais la fraction de français les plus pauvres, qui totalisent 10% de la richesse totale. La bande la plus à droite, ce sont les français les plus riches, qui accumulent 10% de la richesse totale. On peut d’ailleurs être un peu plus éclairant (je pense) en mettant des bandes non pas tous les 10% du revenu cumulé, mais faire l’équivalent du premier graphique, avec 10% des français. Bien entendu, les 10% les plus pauvres possède un très faible pourcentage de la richesse totale. En fait, les 30% les plus pauvres possèdent, tous ensemble, 5% de la richesse totale. Alors que les 10% les plus riches possèdent un tiers de la richesse. C’est ce qu’on peut voir sur le graphique ci-dessous

Ce graphique est – je pense – juste, au sens des normes sémiologiques, et permet de faire passer le message du graphique initial.

Après, on peut tenter une vision qui irait encore plus loin, avec non plus un taux en abscisse, mais un niveau de richesse. Aussi, le taux d’imposition, en fonction du revenu brut mensuel (sur une échelle logarithmique) ressemble au graphique ci-dessous

Forecasting Natural Catastrophes (is rather difficult)

Following my previous post, I wanted to spend more time, on the time series with “global weather-related disaster losses as a proportion of global GDP” over the time period 1990-2016 that Roger Pilke sent me last night.

db=data.frame(year=1990:2016,
ratio=c(.23,.27,.32,.37,.22,.26,.29,.15,.40,.28,.14,.09,.24,.18,.29,.51,.13,.17,.25,.13,.21,.29,.25,.2,.15,.12,.12))

In my previous post, I spend some time explaining that we should provide some sort of ‘confidence interval’ when we try to predict a pattern. That was what we call ‘model uncertainty’. But there are two (important) issues that I did not mention. (1) it is a time series, so why not use techniques dedicated to time series objects ? (2) we do not really care actually about ‘model uncertainty’ (unless we want to assess if a decreasing trend is significant, or not), and we care more about real prediction uncertainty: in the next ten years, what could be the range for the this ratio, with some given probability (say 95%)? Could we say that with 95% chance the global weather-related disaster losses as a proportion of global GDP should be (each year) within 0 and 0.35 or 0 and 0.7?

A first idea might be to use exponential smoothing techniques (without a seasonal component here).

ratio=ts(db$ratio,start=1990,frequency=1)
plot(ratio,xlim=c(1990,2030))
hw=HoltWinters(ratio,gamma=FALSE)
phw=predict(hw,n.ahead=15,prediction.interval = TRUE)
plot(hw,phw,xlim=c(1990,2030))
polygon(c(2017:2031,rev(2017:2031)), c(phw[,2],rev(phw[,3])),border=NA,col=rgb(0,0,1,.2))

The decreasing trend is coming from the fact that exponential smoothing is here a linear regression, with weight exponentially decaying with time (the older, the smaller the weight). But we cannot use that prediction, since the ratio cannot (obviously) be negative. So why not consider, here, the logarithm of the ratio

plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
hw=HoltWinters(log(ratio),gamma=FALSE)
phw=predict(hw,n.ahead=15,prediction.interval = TRUE)
abline(v=2016,lty=2,col="grey")
lines(2017:2031,exp(phw[,2]),col="blue")
lines(2017:2031,exp(phw[,3]),col="blue")
lines(c(1992:2016,2017:2031),c(exp(hw$fitted[,1]),exp(phw[,1])),col="red")
polygon(c(2017:2031,rev(2017:2031)),exp(c(phw[,2],rev(phw[,3]))),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

The confidence band is huge, here. What if we consider some ARIMA model here?

fit=auto.arima(ratio)
farma=forecast(fit,15)
farma=cbind(as.numeric(farma$fitted)[1:15],as.numeric(farma$lower[,1]),as.numeric(farma$upper[,1]),as.numeric(farma$lower[,2]),as.numeric(farma$upper[,2]))
plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
abline(v=2016,lty=2,col="grey")
lines(2017:2031,farma[,4],col="blue")
lines(2017:2031,farma[,5],col="blue")
lines(2017:2031,farma[,1],col="red")
polygon(c(2017:2031,rev(2017:2031)),c(farma[,4],rev(farma[,5])),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

Here, there is an intercept, but no dynamics for the time series (which is considered, here, as a pure white noise). We get exactly the same if we consider the average value of the series

fit=lm(ratio~1,data=db)
s=summary(fit)$sigma
plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
abline(v=2016,lty=2,col="grey")
ndb=data.frame(year=2017:2031)
pf=predict(fit,newdata=ndb)
farma=cbind(pf,pf-1.96*s,pf+1.96*s)
lines(2017:2031,farma[,2],col="blue")
lines(2017:2031,farma[,3],col="blue")
lines(1990:2031,c(predict(fit),farma[,1]),col="red")
polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

Here, we get back to my previous post, if we want to consider a possible trend (and not only an intercept)

fit=lm(ratio~year,data=db)
s=summary(fit)$sigma
plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
abline(v=2016,lty=2,col="grey")
ndb=data.frame(year=2017:2031)
pf=predict(fit,newdata=ndb)
farma=cbind(pf,pf-1.96*s,pf+1.96*s)
lines(2017:2031,farma[,2],col="blue")
lines(2017:2031,farma[,3],col="blue")
lines(1990:2031,c(predict(fit),farma[,1]),col="red")
polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

Again, the confidence region is not based on inference related error, but on model uncertainty: we try to visualize where future observations might be with (say) 95% chance. Note we can also consider (why not?) a quadratic regression

fit=lm(ratio~poly(year,2),data=db)
s=summary(fit)$sigma
plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
abline(v=2016,lty=2,col="grey")
ndb=data.frame(year=2017:2031)
pf=predict(fit,newdata=ndb)
farma=cbind(pf,pf-1.96*s,pf+1.96*s)
lines(2017:2031,farma[,2],col="blue")
lines(2017:2031,farma[,3],col="blue")
lines(1990:2031,c(predict(fit),farma[,1]),col="red")
polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

I am usually not a huge fan of those polynomial regression, but recently, I’ve seen that a lot in economic papers (like “if it’s not linear, add a squared version of the explanatory variable”, which is a rather odd strategy, I’ll publish some posts on that issue this year).

Here again, it might be more clever to consider a logarithmic transformation of the ratio, to insure that the ratio remains positive

fit=lm(log(ratio)~year,data=db)
s=summary(fit)$sigma
plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
abline(v=2016,lty=2,col="grey")
ndb=data.frame(year=2017:2031)
pf=predict(fit,newdata=ndb)
farma=cbind(exp(pf+s^2/2),exp(pf-1.96*s),exp(pf+1.96*s))
lines(2017:2031,farma[,2],col="blue")
lines(2017:2031,farma[,3],col="blue")
lines(1990:2031,c(exp(predict(fit)+s^2/2),farma[,1]),col="red")
polygon(c(2017:2031,rev(2017:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

Observe that future trend is mainly driven by the three latest observations, that were rather low (compared with older observations). What if we remove them?

dbna=db
db$ratio[25:27]=NA
fit=lm(ratio~1,data=dbna)
s=summary(fit)$sigma
plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
abline(v=2016-3,lty=2,col="grey")
ndb=data.frame(year=2014:2031)
pf=predict(fit,newdata=ndb)
farma=cbind(pf,pf-1.96*s,pf+1.96*s)
lines(2014:2031,farma[,2],col="blue")
lines(2014:2031,farma[,3],col="blue")
lines(1990:2031,c(predict(fit)[1:24],farma[,1]),col="red")
polygon(c(2014:2031,rev(2014:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

More funny, if we consider a quadratic regression, we obtain an increasing trend for the future

fit=lm(ratio~poly(year,2),data=dbna)
s=summary(fit)$sigma
plot(db$year,db$ratio,type="l",xlim=c(1990,2030),ylim=c(-.2,.7),xlab="year",ylab="ratio")
abline(v=2016-3,lty=2,col="grey")
ndb=data.frame(year=2014:2031)
pf=predict(fit,newdata=ndb)
farma=cbind(pf,pf-1.96*s,pf+1.96*s)
lines(2014:2031,farma[,2],col="blue")
lines(2014:2031,farma[,3],col="blue")
lines(1990:2031,c(predict(fit)[1:24],farma[,1]),col="red")
polygon(c(2014:2031,rev(2014:2031)),c(farma[,2],rev(farma[,3])),border=NA,col=rgb(0,0,1,.2))
abline(h=0,lty=2)

As we can see, it is rather difficult to get relevant prediction for the future, based on 25 observations…. If anyone has a suggestion, comments are open…

 

What is a Linear Trend, by the way?

I had a very strange discussion on twitter (yes, another one), about regression curves. I think it started with a tweet based on some xkcd picture (just for fun, because it was New Year’s Day)

There were comments on that picture, by econometricians, mainly about ‘significant’ trends when datasets are very noisy. And I mentioned a graph that I saw earlier, a couple of days ago

Let us reproduce that graph (Roger kindly sent me the dataset)

db=data.frame(year=1990:2016,
ratio=c(.23,.27,.32,.37,.22,.26,.29,.15,.40,.28,.14,.09,.24,.18,.29,.51,.13,.17,.25,.13,.21,.29,.25,.2,.15,.12,.12))
library(ggplot2)

The graph is here (with the same aesthetic conventions as Roger’s initial graph, i.e. using some sort of barplot)

ggplot(db, aes(year, ratio)) +
geom_bar(stat="identity") +
stat_smooth(method = "lm", se = FALSE)

My point was that we miss the ‘confidence band’ of the regression

In R, at least, it is quite natural to get (and actually, it is the default version of the graph function)

ggplot(db, aes(year, ratio)) +
geom_bar(stat="identity") +
stat_smooth(method = "lm", se = TRUE)

It is hard to claim that the ‘regression line’ is significant (in the sense “significantly non horizontal”). To be more specific, if we look at the output of the regression model, we get

summary(lm(ratio~year,data=db))

Coefficients:
Estimate    Std. Error t value Pr(>|t|)
(Intercept) 9.158531 4.549672  2.013 0.055 .
year       -0.004457 0.002271 -1.962 0.061 .
---
Signif. codes: 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(which is exactly what Roger used in his graph to plot his red straight line). The p-value of the estimator of the slope, in a linear regression model is here 6%. But I found Roger’s point puzzling

See also

First of all, let us get back to a more standard graph, with a scatterplot, and not bars,

ggplot(db, aes(year, ratio)) +
stat_smooth(method = "lm") +
geom_point()

Here, we observe points https://latex.codecogs.com/gif.latex?\{y_{1990},y_{1991},\cdots,y_{2016}\}. In order to draw that blue line, we assume (Econometrics 101, actually) that those observations are realizations of random variables https://latex.codecogs.com/gif.latex?\{Y_{1990},Y_{1991},\cdots,Y_{2016}\}. Randomness here does not come from a survey, or from ‘balls in an urn’. Randomness is because hurricanes and floods are themselves seen are realizations of random events. Yes, there might be measurement errors, but that’s not where randomness comes from (here). When we talk about ‘randomness’, it should be related to ‘model error’ i.e. the error we make if we consider a linear model (here), that is

https://latex.codecogs.com/gif.latex?Y_t=\beta_0+\beta_1t+\varepsilon_t

Even if observations are not obtained from balls in an urn, there is some kind of randomness here. Randomness means that we might have errors (random errors) around the estimated value (that is on the blue curve), https://latex.codecogs.com/gif.latex?y_t=\widehat{y}_t+\widehat{\varepsilon}_t. One might consider a nonlinear model to reduce the error,

ggplot(db, aes(year, ratio)) +
geom_point() +
geom_smooth()

but in the case, the danger is to overfit

So yes, when we fit a linear model, there is always some kind of randomness, and it is possible to get a ‘confidence band’, that will be very useful for predictions (e.g. for reinsurance purpose here).

Votes Populaires et Grands Électeurs

La semaine dernière, j’avais mis un billet en ligne où je m’interrogeais sur la pertinence de faire commencer l’axe des ordonnées en 0, quand on regardait les élections américaines, en représentant le nombre de votes obtenus. Mon point était que 0 n’avait pas vraiment de sens, quand on regardait républicains et démocrates (en 2016 en tous cas) car ces deux parties obtiendront toujours plusieurs dizaines de millions de voix, quoi qu’il se passe.

Mais les votes populaires importent finalement peu, compte tenu du mode de scrutin. Le plus importants est le nombre de grands électeurs du collège électoral.

> electoral=as.matrix(data.frame(
+ red=c(306,206,173),
+ blue=c(232,332,365)))
> electoral=electoral[3:1,]
> rownames(electoral)=c(2008,2012,2016)
> barplot(t(electoral), col=c("red","darkblue")

Faire commencer le graphique en 0 a ici davantage de sens, car il serait possible d’imaginer qu’un candidat obtienne très peu de grands électeurs. Et force est de constater que ce graphique est  parlant. On peut même tenter d’aller un peu plus loin, en superposant les deux… Le premier soucis est qu’on ne peut pas superposer deux graphiques en barres, donc il faut revenir aux lignes (ou aux points, mais comme c’est une évolution dans le temps, on peut les relier). Le second soucis est que les deux n’ont pas la même échelle (pour les grands électeurs, on en a en tout 539, alors que la population américaine se compte en dizaine de millions d’habitants). On pourrait par exemple tracer

ou encore (de manière toute aussi valide)

Bref, les doubles axes posent soucis (mais ce n’est pas nouveau). La solution serait peut être de travailler en pourcentage. Pour les grands électeurs, on peut mettre le pourcentage obtenu parmi les 539 électeurs. Pour les nombres de votes populaires, on peut se contenter (je pense) de comparer les deux (i.e. pourcentage de voix obtenues par le rouge par rapport au nombre voix obtenues par le rouge et le bleu)

C’est plus parlant, non ?

Les gens bien élevés font commencer les barres en zéro

Il y a quelques jours, je mentionnais sur twitter un graphique mis en ligne sur imgur.com, au sujet des trois dernières élections présidentielles américaines,

Le premier point est que le graphique a été fait alors que les dépouillements n’étaient pas encore finis. 72 heures après la clôture des bureaux de votes, plusieurs millions de bulletins n’avaient pas été dépouillés, comme me l’ont fait noté François Gourio et Tom Roud. Comme on visualise ici le nombre de votants, et pas le pourcentage de votants, c’est génant… Le second point avait été souligné par Marin Dacos, et par plusieurs personnes, toujours sur twitter, qui disaient (je résume un peu)

On ne peut pas faire des diagrammes en barres qui ne commencent pas en zéro ! Sacrilège !

C’est en effet une règle qui semble bien établie. Enfin, quand je dis ‘règle‘, on pourrait penser qu’il s’agit d’une ‘suggestion‘. Des sortes de ‘règles de savoir vivre’, à la Nadine de Rothschild,

Quand on est bien élevé, il est de bon ton de faire commencer ses diagrammes en barres en zéro.

Et cette règle semble importante, car c’est la première que donne Nathan Yau dans ses règles pour faire des graphiques

A bar chart’s baseline must start at zero.

L’intuition est que dans un diagramme en barres, la longueur des barres a un sens.

The bar chart relies on length to show data. Shorter bars represent lower values, and longer bars represent greater values. Compare bar lengths to compare values. That’s how it works. When you shift the baseline, you distort the visual.

Quittons nos élections un instant. Regardons par exemple les temps mis au 100 mètres lors des Championnats du Monde d’Athlétisme.

> year=c(1983,1987,seq(1991,2015,by=2))
> time=c(10.07,9.93,9.86,9.87,9.97,9.86,9.80,9.82, 10.07,9.88,9.85,9.58,9.92,9.77,9.79)
> names(time)=year
> barplot(time,col="darkgreen",beside=TRUE)

Oui, on a compris, les barres ont une interpréation physique, et il ne faut pas les raccourcir ! C’est le temps mis par un sportif pour courrir 100 mètres. Mais force est de constater qu’on ne voit pas grand chose… Allez, souvenez vous en 2009,

Sur notre diagramme en barre, on ne voit rien d’exceptionnel… Allez, osons faire commencer l’axe des ordonnées à 9 secondes,

> barplot(time,ylim=c(9,10.2),col="darkgreen")

Ok, je suis un hérétique… Certes, je mesure ici le temps au delà de 9 secondes, ce qui ne veut pas dire grand chose… J’opère ici la “distortion” dont parle Nathan Yau…. mais j’ai quand même l’impression qu’on voit mieux l’exploit de Bolt (sans laquelle on ne voit absolument rien). Cela dit, on pourrait même tenter quelque chose d’encore plus violent

> barplot(time,ylim=c(min(time),10.2), col="darkgreen")

C’est ridicule comme graphique n’est-ce pas ? Cela dit, à peu de chose près, ce graphique peut se lire en translatant, et en disant qu’il s’agit le temps mis en plus, par rapport au record de Bolt en 2009,

et ce graphique, lui, a du sens et suit la règle de Nathan Yau : il commence en 0.

Sours R, naturellement, un diagramme en baton commencera en 0  mais ça ne sera pas le cas pour une courbe,

Ici, le graphique commence à peu près au minimum de nos observations. Et on comprend avec ce graphique le caractère exceptionnel du temps mis par Bolt, en 2009. Faire commencer le graphique en 0 est ridicule,

C’est presque pire lorsqu’on représente l’évolution du record du monde, au 100 mètres,

Ca n’a pas d’allure ! On ne voit rien sur ce graphique. Alors qu’en changeant l’axe des ordonnées, ça devient plus clair !

Bref, je suis surpris qu’on passe pour un hérétique quand on dit que ce n’est pas génant de ne pas commencer en 0. Ca dépend du message que l’on souhaite faire passer avec le graphique.

En fait, j’ose affirmer que commencer en 0 n’a pas sens ici, car il est impossible qu’un être humain court le 100 mètres en 0 secondes. Cette référence n’a pas lieu d’être ! Mais supposons qu’on nous dise que pour des raisons morphologiques, on pourrait imaginer qu’un homme court, un jour, le 100 mètres en 9 secondes, mais qu’il serait physiologiquement impossible de passer en dessous. Dans ce cas, j’ai envie de croire que la réprésentation en commençant à 9 secondes aurait du sens…

On voit ici non seulement la différence entre Bolt en 2009 et les autres champions du monde, mais on voit aussi la distance qui sépare ces athlètes d’un record absolu. Je suis désolé, mais je ne trouve pas ce graphique dérangeant. Certes, il convient de lire (rapidement) l’axe des ordonnées, pour éviter d’affirmer que Bolt aurait mis 2 fois moins de temps en 2009 que le champion du monde de 2003 pour courrir le 100 mètres (car la barre semble deux fois plus haute).

Revenons maintenant à nos élections américaines, en reprenant les nombres de votes obtenus par candidats en 2008, 2012 et 2016 (avec des chiffres mis à jour)

> votes=data.frame(red=c(59948323,60933504, 62213384),blue=c(69498516,65915795,64154610))
> MV=as.matrix(t(votes))
> colnames(MV)=c(2008,2012,2016)
> barplot(MV/1e6, main="",
+         xlab="", col=c("red","darkblue"),
+         beside=TRUE)

Là encore, j’aurais presque envie d’affirmer que commencer en 0 n’a pas grand sens. Tout d’abord parce qu’il est dur d’imaginer que le candidat des Républicains comme le candidat des Démocrates ait un score proche de zéro. Ensuite, il est aussi difficile d’imaginer qu’un candidat ait deux fois plus de vote que l’autre. Car on peut comparer des barres quand une des deux fois plus grande que l’autre, pas quand elle est 14% plus grande… ça ne se voit pas à l’oeil nu… disons que la différence entre 14% et 18% n’est pas visible. Aussi, quand on connait le contexte, on sait qu’un candidat ne peut pas avoir deux fois plus que l’autre… Donc quand on voit un diagramme à baton de ce genre, on va regarder les valeurs des ordonnées.

> barplot(MV/1e6, main="",
+         xlab="", col=c("red","darkblue"),
+         beside=TRUE,ylim=c(52,72))

J’ai envie de dire qu’on pourrait accepter de fausser la représentation, parce que ce n’est pas ce qu’on veut illustrer avec ce graphique. Ce qu’on veut montrer c’est que le nombre de votes rouges varie augmente un peu, et que le nombre de votes bleus a baissé. Avec des courbes (là on a le droit), on aurait la représentation suivantes

mais je trouve le diagramme en barres plus joli. Et je trouve dommage de dire que le diagramme avec les barres juste au dessus est faussé, alors que le diagramme avec les courbes ne l’est pas… Si on fait commencer en zéro, comme le dit la règle, on ne voit plus rien. Sauf qui ça se joue dans un mouchoir de poche,

Là encore, si on veut respecter la règle du zero, et faire malgré tout des barres, peut être faut-il représenter les variations du nombre de votants (en millions), d’une élection sur l’autre

mais je n’aime pas trop…. Là encore, le graphique le plus parlant reste celui mentionné dans mon tweet

(la nuance est peut être que je ne mets pas ici de barre horizontale en bas). On peut voir le message qu’on veut faire passer avec le graphique de manière beaucoup plus claire. Certains diraient que je fausse l’interprétation, mais je trouve l’affirmation excessive. Il s’agit plutôt de déformer pour s’assurer justement qu’un message (juste et vrai) passe mieux. C’est ce qu’on voyait avec l’exemple du 100 mètres : sans déformation, le seul message que l’on retient est que le record de Bolt est anecdotique.

J’aurais presque envie de croire que la règle ne s’applique pas si la valeur zéro n’est pas atteignable (et n’aurait pas de sens). Mais je ne demande pas mieux que d’être convaincu que même ici, la règle de commencer les barres en zéro doit s’appliquer !

Dominance stochastique, ordre 1 ou ordre 2 ?

Mercredi dernier, on commençait le cours avec une présentation rapide du test de Wilcoxon-Mann-Whitney (parfois appelé test de la somme des rangs, de Wilcoxon). Ce test a été proposé par Frank Wilcoxon en 1945, et (presque en même temps) par Henry Mann et Donald Whitney (en 1947), dans un article intitulé On a Test of Whether one of Two Random Variables is Stochastically Larger than the Other. Le titre est clair. Ce test (non paramétrique) a vocation à tester l’hypothèse selon laquelle la distribution dans deux échantillons est la même. Ou, pour reprendre les notations du derniers billets, la distribution est la même dans deux groupes.

Autrement dit, on va chercher à tester, ici, https://latex.codecogs.com/gif.latex?H_0:F=G, au sens où https://latex.codecogs.com/gif.latex?H_0:\forall%20x,F(x)=G(x). Dans le test de Kolmogorov Smirnov, l’hypothèse alternative est https://latex.codecogs.com/gif.latex?H_1:\exists%20x,F(x)\neq%20G(x).

Ici, comme le notent d’ailleur Mann et Whitney, “the test is shown to be consistent with respect to the class of alternatives https://latex.codecogs.com/gif.latex?F(x)%3EG(x) for every https://latex.codecogs.com/gif.latex?x“. Cette propriété est un peu particulière, et correspond à la notion de dominance stochastique à l’ordre 1. https://latex.codecogs.com/gif.latex?F\preceq_1%20Gcorrespond à https://latex.codecogs.com/gif.latex?F(x)\geq%20G(x) https://latex.codecogs.com/gif.latex?\forall%20x, ou encore, pour tout fonction https://latex.codecogs.com/gif.latex?u(\cdot) croissante, si https://latex.codecogs.com/gif.latex?\mathbb{E}[u(X)]\leq\mathbb{E}[u(Y)] avec https://latex.codecogs.com/gif.latex?X\sim%20F et https://latex.codecogs.com/gif.latex?Y\sim%20G.

Mais regardons un peu le test. On considère deux populations https://latex.codecogs.com/gif.latex?\boldsymbol{x}=\{x_1,\cdots,x_{n_x}\} et https://latex.codecogs.com/gif.latex?\boldsymbol{y}=\{y_1,\cdots,y_{n_y}\}. L’idée est de réordonner les https://latex.codecogs.com/gif.latex?(n_x%20+%20n_y) éléments de https://latex.codecogs.com/gif.latex?\boldsymbol{z}=\boldsymbol{x}\cup\boldsymbol{y} par ordre croissant. On note https://latex.codecogs.com/gif.latex?\boldsymbol{r} le vecteur des rangs associés. Soit https://latex.codecogs.com/gif.latex?S_x la somme des https://latex.codecogs.com/gif.latex?n_x (premiers) rangs des éléments de https://latex.codecogs.com/gif.latex?\boldsymbol{x}, i.e. https://latex.codecogs.com/gif.latex?S_x=r_1+\cdots+r_{n_x}.

Comme l’ont montré Wilcoxon, Mann et Whitney, sous https://latex.codecogs.com/gif.latex?H_0, https://latex.codecogs.com/gif.latex?S_x suit une “loi connue”, où classiquement, “connue” veut dire tabulée pour de petits échantillons, et pour des échantillons plus grands.

Histoire d’illustrer, considerons deux échantillons, tirés suivant deux lois https://latex.codecogs.com/gif.latex?F et https://latex.codecogs.com/gif.latex?G avec https://latex.codecogs.com/gif.latex?F\preceq_1%20G, par exemple deux lois normales de moyennes différentes (et de même variance)

> u=seq(-3,4,length=251)
> plot(u,pnorm(u,0,1),col="red",type="l")
> lines(u,pnorm(u,.5,1),col="darkblue")

Effectivement, nos deux fonctions de répartition sont ordonnées, et on est bien dans une situation de dominance stochastique à l’ordre 1. Les densité sont ici

>plot(u,dnorm(u,0,1),col="red",type="l")
>lines(u,dnorm(u,.5,1),col="darkblue")

On peut générer 10,000 paires d’échantillons, de taille 100 (respectivement), et regarder la distribution de https://latex.codecogs.com/gif.latex?S_x (et celle de https://latex.codecogs.com/gif.latex?S_y).

> P=matrix(NA,1e4,3)
> for(i in 1:1e4){
+  x=c(rnorm(100,0,1),rnorm(100,.5,1))
+  P[i,1:2]=c(sum(rank(x)[1:100]),sum(rank(x)[101:200]))
+  x0=c(rnorm(100,0,1),rnorm(100,0,1))
+  P[i,3]=  sum(rank(x0)[1:100])
+}
> plot(density(P[,1]),col="red"))
> lines(density(P[,2]),col="darkblue")
> lines(density(P[,3]),lty=2)

Sous https://latex.codecogs.com/gif.latex?H_0, la somme aurait eu la distribution avec les pointillés. On en est manifestement assez loin, autrement dit, on va rejeter ici l’hypothèse de lois égales, https://latex.codecogs.com/gif.latex?H_0:F=G.

Intuitivement, c’est normal, car effectivement, les rangs dans un cas sont – en moyenne – plus faibles que pour l’autre,

Au lieu de simuler des lois dominées, à l’ordre 1, on peut tenter des lois dominées à l’ordre 2. Rappelons que https://latex.codecogs.com/gif.latex?F\preceq_2%20G signifie

https://latex.codecogs.com/gif.latex?\int_{-\infty}^x%20F(t)dt%20\geq\int_{-\infty}^x%20G(t)dt,\forall%20x

ou encorehttps://latex.codecogs.com/gif.latex?\mathbb{E}[u(X)]\leq\mathbb{E}[u(Y)]

avec https://latex.codecogs.com/gif.latex?X\sim%20F et https://latex.codecogs.com/gif.latex?Y\sim%20G. où https://latex.codecogs.com/gif.latex?u(\cdot) est unne fonction croissante et concave.

En gestion des risques (ou sur les inégalités), on va supposer en plus que https://latex.codecogs.com/gif.latex?\mathbb{E}[X]=\mathbb{E}[Y], ce qui correspond à une notion de “mean preserving spread” (on peut relire le document sur la notion de
single crossing, par Alain Chateauneuf). Une caractérisation est

https://latex.codecogs.com/gif.latex?\mathbb{E}[u(X)]\leq\mathbb{E}[u(Y)]https://latex.codecogs.com/gif.latex?X\sim%20F et https://latex.codecogs.com/gif.latex?Y\sim%20G. où https://latex.codecogs.com/gif.latex?u(\cdot) est concave. Une caractérisation alternative est https://latex.codecogs.com/gif.latex?Y=X+\varepsilon avec https://latex.codecogs.com/gif.latex?\mathbb{E}[\varepsilon\vert%20X]=0 (on pourra relire Rothschild et Stiglitz, Increasing Risk: A Definition)

Dans le cas des lois normales, on a cette notion de ‘mean preserving spread’ lorsque les moyennes sont identiques, mais pas les variances. On a alors plus de dispersion,

> u=seq(-4,4,length=251)
> plot(u,pnorm(u,0,1),col="red",type="l")
> lines(u,pnorm(u,0,1.5),col="darkblue")

Effectivement, nos deux fonctions de répartition sont ordonnées, et on est bien dans une situation de dominance stochastique à l’ordre 2, avec du single crossing. Les densités sont ici

> plot(u,dnorm(u,0,1),col="red",type="l")
> lines(u,dnorm(u,0,1.5),col="darkblue")

On peut générer 10,000 paires d’échantillons, de taille 100 (respectivement), et regarder la distribution de https://latex.codecogs.com/gif.latex?S_x (et celle de https://latex.codecogs.com/gif.latex?S_y).

> P=matrix(NA,1e4,3)
> for(i in 1:1e4){
+  x=c(rnorm(100,0,1),rnorm(100,0,1.5))
+  P[i,1:2]=c(sum(rank(x)[1:100]), sum(rank(x)[101:200]))
+  x0=c(rnorm(100,0,1),rnorm(100,0,1))
+  P[i,3]=  sum(rank(x0)[1:100])
+}
> plot(density(P[,1]),col="red"))
> lines(density(P[,2]),col="darkblue")
> lines(density(P[,3]),lty=2)

Autrement dit, ici les sommes des rangs sont identiques, avec les deux échantillons,  et donc on aura tendance à accepter ici l’hypothèse de lois égales, https://latex.codecogs.com/gif.latex?H_0:F=G. L’idée naturelle est que avec les dans le premiers cas, effectivement, on avait plus de petits rangs avec le petit échantillon. Mais plus ici: le fait d’avoir plus de dispersion fait qu’on a en même temps plus de petits, et plus de grands. Et en moyenne, ça ne change rien

Modéliser et prévoir, entre ‘avoir faux’ et ‘avoir pas de chance’

Depuis plusieurs heures, il y a une phrase qui revient sans cesse, et qui m’agace au plus haut point, sur “le modèle est faux…“. C’est une phrase que j’ai entendu pas plus tard qu’au séminaire ce midi, lorsqu’un collègue a dit que très clairement les modèles étaient faux puisqu’ils n’avaient pas prédit le vainqueur des élections américaines. Par exemple, le jour des élections, Huffington Post annonçait que Donald Trump avait 1,6% de chances d’être élu…

et pourtant, il a gagné. C’est donc pourtant que le modèle était “faux”, c’est évident, non ?

Pourtant, tout ceux qui ont fait un minimum de statistique dans leurs vies le savent… les évènements rares, ça peut se produire,

> sample(c("clinton","trump"),size=1,
+ prob=c(.981,.016))
[1] "clinton"
> sample(c("clinton","trump"),size=1,
+ prob=c(.981,.016))
[1] "clinton"
> sample(c("clinton","trump"),size=1,
+ prob=c(.981,.016))
[1] "trump"

Ce type de raisonnement me dérange car il dénote une mauvaise compréhense de ce qu’est une probabilité. Autre exemple: supposons que j’annonce que ma belle mère a une chance sur un million de mourrir lorsqu’elle prendra la route pendant les prochaines vacances, et qu’elle meurt effectivement cet hiver d’un accident de voiture. Comme pour les élections américaines, on sera surpris, sous le choc. Mais est-ce que ça veut dire que mon modèle était faux ? Non…Si j’avais annoncé qu’elle avait 85% de ne pas passer les vacances en prenant la voiture, j’aurais eu un modèle des plus douteux… Le problème ici c’est qu’on avance des probabilités sur des évènements qu’on ne pourra observer qu’une fois. On ne peut pas le tester suffisement pour savoir s’il est juste, en moyenne.

Cela dit, c’est un problème assez profond, car il est très dur de savoir ce qui relève effectivement de la malchance… C’est en tous cas comme ça que je comprends le strip ci-dessous

Et je trouve dommage que sous prétexte qu’un évènement rare se réalise, on blâme la personne en charge de la prévision… Il y a manifestement une incompréhension profonde de ce que peut être le hasard. C’est triste…

Ce que la courbe ROC (et l’AUC) ne raconte pas

En préparant une intervention pour mardi prochain, j’épluchais les résultats renvoyés pour un exercice, et j’ai eu un résultat assez étrange avec un modèle de classification. J’avais donné la même base cet automne à l’ensae, et j’avais donc près d’une trentaine d’autres modèles, pour comparer (disons plutôt que sur la même base de test, j’ai une trentaine de prévisions). Les observations noires sont celles obtenues cet automne (le trait correspond aux meilleurs AUC sur la base de test), et les observations rouges sont celles que j’ai obtenu pour l’intervention de mardi (là encore, le trait vertical correspond aux meilleurs modèles, au sens du AUC), sur une observation de la base de test,

Ce sont les probabilités prédites (mais j’ai enlevé l’échelle).

Pour presque toutes mes observations, les poids rouges sont bien au dessus des autres… Mais cela n’enlève en rien le fait que l’AUC obtenu (pour les deux modèles rouges) est très bon. C’est effectivement un résultat important (et connu) : le critère AUC (mais plus généralement la courbe ROC en entier) n’indique en rien si la valeur prédite est bonne, ou pas. Il nous dit juste si l’ordre obtenu est correct. Si les valeurs les plus importantes sont effectivement les valeurs pour lesquelles on a un 1, l’AUC sera très bon.

C’est ce qu’on peut observer sur le petit exemple ci-dessous. Considérons un modèle logistique simulé assez simple,

> n=1e3
> set.seed(1)
> x1=rnorm(n)
> x2=runif(n)
> u=-3+x2+x1
> p=exp(u)/(1+exp(u))
> y=rbinom(n,prob=p,size=1)
> library(ROCR)
> df=data.frame(y,x1,x2)
> mean(df$y)
[1] 0.116
> reg=glm(y~.,data=df,family=binomial)
> p=predict(reg,type="response")
> mean(p)
[1] 0.116
> pred1=prediction(p, df$y)
> L=performance(pred1, "tpr", "fpr")

L’AUC est ici

> auc=performance(pred1, "auc")@y.values[[1]]
> auc
[1] 0.7681191

et la courbe ROC est la suivante,

> plot(unlist(L@x.values),unlist(L@y.values),
+ type="s",col="blue")

Supposons que l’on change la constante du modèle logistique,

> reg$coefficients[1]=0

Dans ce cas, notre prévision est assez mauvaise, car la probabilité moyenne prédite est ici

> u=reg$coefficients[1]+reg$coefficients[2]*
+ df$x1+reg$coefficients[3]*df$x2
> p=exp(u)/(1+exp(u))
> mean(p)
[1] 0.6060676

(on est loin des 11,6% de 1 dans la base). Pourtant le AUC est bon

> pred1=prediction(p, df$y)
> L=performance(pred1, "tpr", "fpr")
> auc=performance(pred1, "auc")@y.values[[1]]
> auc
[1] 0.7681191

(c’est en fait la même valeur qu’auparavant, ce qui a du sens puisque la courbe ROC est identique)

> lines(unlist(L@x.values),unlist(L@y.values),
+ type="s",col="red")

Autrement dit, ces outils, classiquement utilisés pour juger la qualité d’un classifieur, ne permettent en aucun cas de dire que la probabilité prédite a du sens. Ces critères permettent juste de dire qu’on identifie assez bien les personnes qui ont le plus de chance d’avoir la réponse 1. Ce qui n’est pas si mal… mais c’est un autre problème que celui d’avoir une probabilité qui soit pertinente.

Estimer une Proportion sur une Question Sensible, par Enquête

Comme le rapportait Slate hier (dans un article qui, malheureusement ne mérite pas d’être lu), un joli sujet a été proposé au baccalauréat au lycée français de Pondichéry, en mathématiques.

Ce sujet est d’autant plus intéressant qu’il revient sur une méthode assez classique pour questionner les gens sur des questions sensibles (ici sur le téléchargement pirate, mais on peut poser la même question sur la fraude, ou des pratiques sexuelles).

Comme le notait , cette technique avait été introduite initialement par Stanley Warner, voilà un demi-siècle,

Cette technique est particulièrement intéressante, tant d’un point de vu intellectuel que pratique. Plusieurs sur Twitter me faisaient remarquer que, contrairement à ce qui était raconté dans Slate, cette technique est effectivement utilisée.

Quand j’enseignais cette technique, je prenais la version simple, avec une pièce de monnaie. Supposons qu’une proportion  de personnes fraudent le fisc. Et qu’une proportion  aime Prince, va à la piscine une fois par mois, a un diplôme de l’enseignement supérieur, etc. Peut importe, la seconde proportion doit pouvoir être estimée sans problème avec une autre enquête (voire par le recensement). Bref, on a un estimateur de . Ici, les personnes interrogées tirent une pièce et répondent soit à la question “avez-vous déjà fraudé au Fisc?” s’ils tombent sur face, et à la question “aimez-vous Prince?” s’ils tombent sur pile. On pourrait dire qu’ils tirent une pièce dans leur tête, mais il faut s’assurer que la probabilité de tomber sur face est effectivement 1/2.

La magie est de noter que la proportion de personnes qui fraudent le fisc peut s’écrire

On peut avoir, par une autre enquête un estimateur du second terme, . Et en regardant des personnes qui ont répondu “oui” à une question inconnue, on a un estimateur pour

Et par différence, on en déduit un estimateur de . Les maths sont dans l’article de Stanley Warner, mais il faut reconnaître que la méthode est vraiment jolie…. Dans le cas du sujet du bac, la mise en oeuvre est un peu différente, mais l’idée est un peu la même….

(si certains ont des exemples pratiques d’application, les commentaires sont ouverts, ou sinon par mail pour joindre des documents en pièce jointe)

Radial Graphs for Time Series

On How to: Weather Radials, there was a nice visualisation of temperatures. Since I am too old fashioned for ggplot2, I wanted to reproduce a similar graph with the old plot style.

Assume that daily temperature is in a vector X (e.g. temperature in Montréal, QC, in 2009). To get a radial plot, use

> n=length(X)
> theta=seq(0,1-1/n,length=n)*2*pi
> r=30+X
> plot(r*cos(pi/2-theta),r*sin(pi/2-theta),type="l",xlab="",ylab="",axes=FALSE)
> for(t in 1:n){
+   if(X[t]>0) CL=rgb(0,0,1,.4)
+   if(X[t]<0) CL=rgb(1,0,0,.4)
+   if(X[t]==0) CL="white"
+   segments((30+X[t])*cos(pi/2-theta[t]),(30+X[t])*sin(pi/2-theta[t]),30*cos(pi/2-theta[t]),30*sin(pi/2-theta[t]),col=CL)
+ }
> for(r in 10*seq(0,6)) lines(r*cos(pi/2-theta),r*sin(pi/2-theta),type="l",col="light blue")

Bonferroni, les anniversaires et les terroristes (ou pas)

Le paradoxe des anniversaires est un exercice de probabilité que l’on découvre assez vite. Mais assez souvent, on l’oublie, alors qu’il est en réalité assez profond. Le problème est le suivant: quel est le nombre de personnes qu’il faut réunir pour avoir (au moins) une chance sur 2 pour que deux personnes aient leur anniversaire le même jour. Si on calcule la probabilité que deux personnes dans un groupe de n fêtent leur anniversaire le même jour, on peut calculer la probabilité que personne ne soit né le même jour, parmi n personnes. Un calcul rapide montre qu’avec 23 personnes, la probabilité qu’au moins deux personnes fête son anniversaire le même jour est 50.73%. Autrement dit, il suffit d’une groupe de 23 personnes (c’était la taille d’une classe d’élèves, du moins à mon époque…).

L’intuition, pour comprendre pourquoi nous avons l’impression que ce nombre est petit est que, naturellement, si on a un événement qui survient avec probabilité p, et qu’on répète l’expérience n fois, la probabilité que rien ne se produise est (1-p)n. Aussi, en première approximation (et j’insiste sur ce point parce que c’est avec ce genre d’argument fallacieux qu’on conclue que la probabilité d’avoir une catastrophe nucléaire dépasse 100%, conduisant à du “certitude statistique“), on en déduit que la probabilité qu’un événement survienne est de l’ordre de np. Sauf qu’il, on compte des paires. Et avec n personnes, on constitue environ n2/2 paires. Et donc la probabilité devient proche de n2p/2… ce qui est bien plus grand que np. C’est grossièrement ce qui explique que notre intuition se sent trompée : avec seulement 23 personnes, on peut constituer plus de 250 paires. La probabilité que toutes les paires correspondent à des dates d’anniversaire différente n’est pas si petit qu’il y parait.

Jeff Ullman parle ainsi de principe de Bonferroni,

“if you look in more places for interesting patterns than your amount of data will support, you are bound to find crap”

Et cette intuition faussée quand on travaille sur des paires se retrouve dans les techniques utilisées par les services de renseigngement, tel que le montre Anand Rajaraman dans Mining Massive DatasetsSupposons qu’il existe un milliard d’individus qui pourraient potentiellement être des terroristes. Que les gens, quels qu’il soient vont à l’hôtel 1 jour sur 100. Un hôtel donné peut accueillir 100 personnes. On suppose qu’il existe 100,000 personnes, de sorte qu’en moyenne ils soient plein tous les jours. Et supposons qu’on regarde sur un historique de 3 ans, i.e. 1000 jours pour avoir des calculs ronds. Supposons que pour traquer les terroristes, on cherche les gens qui ont séjourné, ensemble, dans un hôtel, le même jour, au moins deux fois… On imagine que c’est effectivement rare, et donc suspicieux…

La probabilité que deux personnes données dorment à l’hôtel le même jour est faible, 1/100×1/100, soit 0,0001. Et comme il y a 100,000 hôtel, la probabilité que les deux personnes descendent, en plus, dans le même hôtel, il faut diviser 0,0001 par 100,000. Ce qui donne 10-9 . Mais ça, c’est la probabilité pour qu’ils se croisent une fois ! Si on regarde la probabilité qu’ils se croisent deux fois, en trois ans, on arrive à une probabilité de 10-18 Pour ceux qui ont du mal, ça doit faire environ

0,00000000000000001 %

Et encore. Ici il peut en plus se croiser la seconde fois dans un autre hôtel. La probabilité qu’il se croisent dans le même hôtel que la première fois devient encore plus petite !

Supposons maintenant qu’un organisme de surveillance avisé décider de lister les paires de personnes qui, sur des paires de jours, se sont retrouvé dans le même hôtel. Ça doit être suffisamment rare pour être suspicieux. Rappelons que quand on considère n individus, le nombre de paires possibles est de l’ordre de n2/2. Autrement dit, avec un million de personnes, le nombre de paires est de l’ordre de 5×1017. Et pareil, le nombre de paires de jours est de l’ordre de 5×105. Aussi, le nombre d’événements suspicieux que l’on observe, sur ces trois années est de l’ordre de

(5×1017 ) x (5×10) x10-18 = 250000

Autrement dit, 1/4 de millions de paires de personnes sont des terroristes potentiels ! Soit 1/2 millions d’individus.

Maintenant, supposons qu’il y ait 10 terroristes. S’ils sont assez stupides pour se retrouver dans les mêmes hôtels, les mêmes jours, ça veut dire que sur les 500000 personnes qui vont être surveillées, 10 sont effectivement des terroristes, et on a 499990 faux positifs… Enfin, “faux positifs”… ils ont forcément quelque chose d’autre à se reprocher, non ?

Dynamique de la Pyramide des Ages

Très joli billet sur blog.revolutionanalytics.com avec un code de @kyle_e_walker permettant, très simplement (moyennant une inscription pour avoir une clé permettant d’utiliser l’API du census) de construire une pyramide des âges dynamiques.

> devtools::install_github('walkerke/idbr')
> library(idbr)
> library(ggplot2)
> library(animation)
> library(dplyr)
> library(ggthemes)
> idb_api_key("mykey1239F2f324zf9GGZgege32R2ii4")

On importe alors les données pour les hommes et les femmes,

> male <- idb1('FR', 2010:2050, sex = 'male') %>%
+   mutate(POP = POP * -1,
+   SEX = 'Male')
> female <- idb1('FR', 2010:2050, sex = 'female') %>% mutate(SEX = 'Female')

et on stocke le tout

> france <- rbind(male, female) %>%
+   mutate(abs_pop = abs(POP))

Ensuite, on crée l’animation,

> saveGIF({
+   
+   for (i in 2010:2050) {
+     
+     title <- as.character(i)
+     
+     year_data <- filter(france, time == i)
+     
+     g1 <- ggplot(year_data, aes(x = AGE, y = POP, fill = SEX, width = 1)) +
+       coord_fixed() + 
+       coord_flip() +
+       annotate('text', x = 98, y = -800000, 
+       label = 'Data: US Census Bureau IDB; idbr R package', size = 3) + 
+       geom_bar(data = subset(year_data, SEX == "Female"), stat = "identity") +
+       geom_bar(data = subset(year_data, SEX == "Male"), stat = "identity") +
+       scale_y_continuous(breaks = seq(-1000000, 1000000, 500000),
+       labels = paste0(as.character(c(seq(1, 0, -0.5), c(0.5, 1))), "m"), 
+       limits = c(min(france$POP), max(france$POP))) +
+       theme_economist(base_size = 14) + 
+       scale_fill_manual(values = c('#ff9896', '#d62728')) + 
+       ggtitle(paste0('Population structure of France, ', title)) + 
+       ylab('Population') + 
+       xlab('Age') + 
+       theme(legend.position = "bottom", legend.title = element_blank()) + 
+       guides(fill = guide_legend(reverse = TRUE))
+     print(g1) 
+   }
+ }, movie.name = 'france_pyramid.gif', interval = 0.1, ani.width = 700, ani.height = 600)

Et le résultat est vraiment joli, non ?