Tag Archives: computer

Finding Waldo, a flag on the moon and multiple choice tests, with R

I have to admit, first, that finding Waldo has been a difficult task. And I did not succeed. Neither could I correctly spot his shirt (because actually, it was what I was looking for). You know, that red-and-white striped shirt. I guess it should have been possible to look for Waldo’s face (assuming that his face does not change) but I still have problems with size factor (and resolution issues too). The problem is not that simple. At thehttp://mlsp2009.conwiz.dk/ conference, a price was offered for writing an algorithm in Matlab. And one can even find Mathematica codes online. But most of the those algorithms are based on the idea that we look for similarities with Waldo’s face, as described in problem 3 on http://www1.cs.columbia.edu/~blake/‘s webpage. You can find papers on that problem, e.g. Friencly & Kwan (2009) (based on statistical techniques, but Waldo is here a pretext to discuss other issues actually), or more recently (but more complex) Garg et al. (2011) on matching people in images of crowds.

What about codes in R ? On http://stackoverflow.com/, some ideas can be found (and thank Robert Hijmans for his help on his package). So let us try here to do something, on our own. Consider the following picture,

With the following code (based on the following file) it is possible to import the picture, and to extract the colors (based on an RGB decomposition),

> library(raster)
> waldo=brick(system.file("DepartmentStoreW.grd",
+ package="raster"))
> waldo
class       : RasterBrick
dimensions  : 768, 1024, 786432, 3 (nrow,ncol,ncell,nlayer)
resolution  : 1, 1  (x, y)
extent      : 0, 1024, 0, 768  (xmin, xmax, ymin, ymax)
coord. ref. : NA
values      : C:\R\win-library\raster\DepartmentStoreW.grd
min values  : 0 0 0
max values  : 255 255 255

My strategy is simple: try to spot areas with white and red stripes (horizontal stripes). Note that here, I ran the code on a Windows machine, the package is not working well on Mac. In order to get a better understanding of what could be done, let us start with something much more simple. Like the picture below, with Waldo (and Waldo only). Here, it is possible to extract the three colors (red, green and blue),

> plot(waldo,useRaster=FALSE)

It is possible to extract the red zones (already on the graph above, since red is a primary color), as well as the white ones (green zones on the graphs means a white region on the picture, on the left)

# white component
white = min(waldo[[1]] , waldo[[2]] , waldo[[3]])>220
focalswhite = focal(white, w=3, fun=mean)

# red component
red = (waldo[[1]]>150)&(max(  waldo[[2]] , waldo[[3]])<90)
focalsred = focal(red, w=3, fun=mean)

i.e. here we have the graphs below, with the white regions, and the red ones,

From those two parts, it has been possible to extract the red-and-white stripes from the picture, i.e. some regions that were red above, and white below (or the reverse),

# striped component
striped = red; n=length(values(striped)); h=5
focalsstriped = focal(striped, w=3, fun=mean)

So here, we can easily spot Waldo, i.e. the guy with the red-white stripes (with two different sets of thresholds for the RGB decomposition)

Let us try somthing slightly more complicated, with a zoom on the large picture of the department store (since, to be honest, I know where Waldo is…).

Here again, we can spot the white part (on the left) and the red one (on the right), with some thresholds for the RGB decomposition

Note that we can try to be (much) more selective, playing with threshold. Here, it is not very convincing: I cannot clearly identify the region where Waldo might be (the two graphs below were obtained playing with thresholds)

And if we look at the overall pictures, it is worst. Here are the white zones, and the red ones,

and again, playing with RGB thresholds, I cannot spot Waldo,

Maybe I was a bit optimistic, or ambitious. Let us try something more simple, like finding a flag on the moon. Consider the picture below on the left, and let us see if we can spot an American flag,

Again, on the left, let us identify white areas, and on the right, red ones

Then as before, let us look for horizontal stripes

Waouh, I did it ! That’s one small step for man, one giant leap for R-coders ! Or least for me… So, why might it be interesting to identify areas on pictures ? I mean, I am not Chloe O’Brian, I don’t have to spot flags in a crowd, neither Waldo, nor some terrorists (that might wear striped shirts). This might be fun if you want to give grades for your exams automatically. Consider the two following scans, the template, and a filled copy,

A first step is to identify regions where we expect to find some “red” part (I assume here that students have to use a red pencil). Let us start to check on the template and the filled form if we can identify red areas,

exam = stack("C:\\Users\\exam-blank.png")
red = (exam[[1]]>150)&(max(  exam[[2]] , exam[[3]])<150)
focalsred = focal(red, w=3, fun=mean)
exam = stack("C:\\Users\\exam-filled.png")
red = (exam[[1]]>150)&(max(  exam[[2]] , exam[[3]])<150)
focalsred = focal(red, w=3, fun=mean)

First, we have to identify areas where students have to fill the blanks. So in the template, identify black boxes, and get the coordinates (here manually)

exam = stack("C:\\Users\\exam-blank.png")
black = max(  exam[[1]] ,exam[[2]] , exam[[3]])<50
focalsblack = focal(black, w=3, fun=mean)

The blue points above are where we look for students’ answers. Then, we have to define the vector of correct answers,

points(CORRECTX, CORRECTY,pch=16,col="red",cex=1.3)

Now, let us get back on red areas in the form filled by the student, identified earlier,

exam = stack("C:\\Users\\exam-filled.png")
red = (exam[[1]]>150)&(max(  exam[[2]] , exam[[3]])<150)
focalsred = focal(red, w=5, fun=mean)

Here, we simply have to compare what the student answered with areas where we expect to find some red in,

points(CORRECTX, CORRECTY,pch=1,

Crosses on the graph on the right below are the answers identified as correct (here 13),

> icorrect=values(red)[(750-CORRECTY)*
+ 610+(CORRECTX)]
> points(CORRECTX[icorrect], CORRECTY[icorrect],
+ pch=4,col="black",cex=1.5,lwd=1.5)
> sum(icorrect)
[1] 13

In the case there are negative points for non-correct answer, we can count how many incorrect answers we had. Here 4.

> iuncorrect=values(red)[(750-UNCORRECTY)*610+
> sum(iuncorrect)
[1] 4

So I have not been able to find Waldo, but I least, that will probably save me hours next time I have to mark exams…

Correlations, dimension, and risk measure

Yesterday, while I was attending the IFM2 conference, at HEC Montreal, I heard a nice talk about credit risk, and a comparison between contagion (or at least default correlation), for corporate and retail companies (in the US). And it was mentioned that default correlation was much lower for retail companies than it could be for corporate risk. In a discussion that followed those slides, it was mentioned that banks in the US should actually have been working more with those small firms, since contagion risk was much lower.

A problem here is that the link between correlation, risk and dimension is rather complicated:

  • corporate means a small number of firms, high correlation (and possible large individual losses)
  • retail means a large number of firms (even perhaps extremely large), lower correlation (and small individual losses)

A simple model for default models is based on the assumption that we deal with an exchangeable portfolio (as in a previous post). With the following code, given an (individual) default probability, a default correlation, and a number of firms, it is possible to calculate the probability to have more than a given number of defaults.

 stop.on.error =  FALSE)$value}

a=m*(1-r)/r ;
 for(i in 0:n){
 return(sum(V[1:(x+1)])) }

It is possible to calculate, for a large range of correlations, the probability to have – at least – 20% of default in the portfolio (in order to compare things that are comparable).

for(i in 1:length(A)){

With 20 firms (corporate) we want to have at least 4 defaults, while with 1000 firms (retail) there should be 200 defaults. As mentioned in the previous post, the relationship between correlation and quantiles of sums is not simple. Hence, it might not be monotone. The dotted line is the probability to have at least 4 defaults when default correlation is 50% (around 10%). The plain line is the probability to have at least 200 defaults, as a function of the correlation,


In that case, with only a correlation of 10% among retail firms, the probability of having 20% defaults is the same as the same probability for corporate, but with 50% correlation… One should remember that in portfolio analysis, the links between correlation, dimension and risk measure is a sensitive issue…

Régression médiane et géométrie

Suite à mon exposé d’hier, Pierre me demandait d’argumenter un point que j’avais évoqué oralement sans le justifier: “ une régression médiane passe forcément par deux points du nuage “. Par exemple,


Essayons de le justifier… de manière un peu heuristique. Commençons un peu au hasard, avec une droite qui passe “entre” les points (on admettra que la régression médiane sépare l’espace en deux, avec la moitié des points en dessous, la moitié au dessus, à un près pour des histoires de parité).


c’est joli, mais on peut sûrement faire mieux. En particulier, on va chercher ici à minimiser la somme des valeurs absolue des erreurs (c’est le principe de la régression médiane). Essayons de translater la courbe, vers le haut, ou vers le bas,

for(i in seq(-2,3,by=.5)){

Si on regarde ce que vaut la somme des valeurs absolues des erreurs, on a

d=function(h) sum(abs(y-(h+3.2*x)))

Formellement, l’optimum est ici

> optimize(d,lower=-5,upper=5)
[1] -0.2000091

[1] 2.200009

Bref, retenons cette courbe, et notons qu’elle passe par un des points,

Et c’est assez normal. On commençait avec deux points au dessus, et un en dessous. Soit la somme des valeurs absolues des erreurs . Si on translate de , on passe de à (tant que l’on a toujours deux points au dessus, car celui en dessous sera toujours en dessous). Si on translate de , (là encore tant qu’il n’y a qu’un point en dessous). Bref, on a intérêt à translater vers le haut. Si on dépasse le premier point rencontré, on se retrouve dans la situation inverse, avec un point au dessus et deux au dessous, et on a intérêt à redescendre. Bref, la translation optimale revient a s’arrêter dès qu’on croise un point. Autrement dit, la régression passe forcément par un point. Au moins.

Maintenant, essayons de faire pivoter la courbe autour de ce point, là encore afin de minimiser la somme des valeurs absolues des erreurs,

for(i in seq(-1,5,by=.25)){

La distance est alors, en fonction de la pente de la droite

d2=function(h) sum(abs(y-((y[1]-h*x[1])+h*x)))

Là encore on peut formaliser un peu

> optimize(d2,lower=-5,upper=5)
[1] 2.500018

[1] 1.500018

Et si on regarde cette dernière droite, on passe par deux points,


Pour comprendre pourquoi, comparons les deux cas,

  • en faisant un pivot vers le bas, pour un des points, la valeur absolue des erreurs augmente alors que pour l’autre, elle diminue

  • en faisant un pivot vers le haut, là encore pour un des points, la valeur absolue des erreurs augmente alors que pour l’autre, elle diminue, mais en sens inverse par rapport au cas précédent,

Et pour faire simple, dans le premier cas, le gain (sur la baisse) compense la perte, alors que dans le second cas, c’est le contraire. Bref, on a intérêt à pivoter vers le bas, ou vers le point à droite. Jusqu’à l’atteindre. Et optimalement, on passera par ce point. Bref, on passera par deux points. Et je laisse les plus courageux regarder avec plus de trois points, mais c’est toujours pareil…

Sondages et prévisions de séries temporelles

Ce 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


(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")
for(k in 2:4){

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


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

for(k in 2:4){

avec, là encore, la commande suivante pour lisser


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),


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,

for(k in c(1,3:4)){

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


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)


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

variance= -proba%*%t(proba)
for(i in 1:nrow(M)){
for(k in 1:4){

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

for(k in 2:4){

(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é),

for(k in 2:4){

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


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,

ETS = ets(X)
for(k in 2:4){

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,

for(sim in 1:20){

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

for(sim in 1:1000){

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 !

Sur le lissage exponentiel

Sur le lissage exponentiel, je pourrais renvoyer vers Gardner (1985), et la version remise au goût du jour, Gardner (2005). J’avais pris comme position, dans le cours, de présenter rapidement le lissage exponentiel, en notant qu’on les évoquerais à nouveau plus tard comme cas particulier des modèles ARIMA. Comme le note la dernière version, Gardner (2005), ce n’est pas forcément l’unique manière de voir, et le cours pourrait être orienté complètement autour des notions de lissage exponentiel (c’est d’ailleurs le point du vue de livre Hyndman, Koehler, Ord & Snyder (2008)), “when Gardner (2005) appeared, many believed that exponential smoothing should be disregarded because it was either a special case of ARIMA modeling or an ad hoc procedure with no statistical rationale. As McKenzie (1985) observed, this opinion was expressed in numerous references to my paper. Since 1985, the special case argument has been turned on its head, and today we know that exponential smoothing methods are optimal for a very general class of state-space models that is in fact broader than the ARIMA class.”

N’étant toujours pas convaincu, j’évoquerais à nouveau le lissage exponentiel quand on présentera les modèles ARIMA. En attendant, un peu de code pour mieux comprendre ce qui est fait quand on fait du lissage exponentiel.

Commençons par le lissage exponentiel simple, i.e.


https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/le04.gif désigne un poids attribué à la nouvelle observation dans la fonction de lissage https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex05.gif. Le code pour lisser une série est alors assez simple,

> library(datasets)
> X=as.numeric(Nile)
> Lissage=function(a){
+  T=length(X)
+  L=rep(NA,T)
+  L[1]=X[1]
+  for(t in 2:T){L[t]=a*X[t]+(1-a)*L[t-1]}
+  return(L)
+ }
> plot(X,type="b",cex=.6)
> lines(Lissage(.2),col="red")


Sur la figure ci-dessus, non visualise l’impact de https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse16.gif sur le lissage. La prévision que l’on fait à la date https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lex14.gif, pour un horizon https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse12.gif est alors https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse11.gif. Il est alors possible de voir le poids comme un paramètre et on va alors chercher le poids optimal. La stratégie classique est de minimiser l’erreur de prédiction commise à un horizon de 1


> V=function(a){
+  T=length(X)
+  L=erreur=rep(NA,T)
+  erreur[1]=0
+  L[1]=X[1]
+  for(t in 2:T){
+  L[t]=a*X[t]+(1-a)*L[t-1]
+  erreur[t]=X[t]-L[t-1] }
+ return(sum(erreur^2))
+ }

Ici, on obtient comme poids optimal

> A=seq(0,1,by=.02)
> Ax=Vectorize(V)(A)
> plot(A,Ax,ylim=c(min(Ax),min(Ax)*1.05))
> optimize(V,c(0,.5))$minimum
[1] 0.246581

On notera que c’est ce que suggère la fonction de R,

> hw=HoltWinters(X,beta=FALSE,gamma=FALSE,l.start=X[1])
> hw
Holt-Winters exponential smoothing without trend an seasonal comp.

HoltWinters(x = X, beta = FALSE, gamma = FALSE, l.start = X[1])

Smoothing parameters:
alpha:  0.2465579
beta :  FALSE
gamma:  FALSE

a 805.0389

> plot(hw)
> points(2:(length(X)+1),Vectorize(Lissage)(.2465),col="blue")

Dans le cas du lissage exponentiel double, i.e.




Dans ce cas, la prédiction est https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse13.gif. Le code pour faire un lissage double est là encore assez simple,

> Lissage=function(a,b){
+  T=length(X)
+  L=B=rep(NA,T)
+  L[1]=X[1]; B[1]=0
+  for(t in 2:T){
+  L[t]=a*X[t]+(1-a)*(L[t-1]+B[t-1])
+  B[t]=b*(L[t]-L[t-1])+(1-b)*B[t-1] }
+ return(L)
+ }

Sur la figure suivante, on visualise l’évolution du lissage en fonction de https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse15.gif (https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse16.gif étant ici fixé),
(le lissage simple – avec le même poids https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lisse16.gif – apparaissant ici en trait clair).

Visualisation des résidus pour des données spatiales

Mardi, nous allons travailler un peu sur la modélisation du nombre d’homicides aux États-Unis, à partir de la base

> US=read.table("http://freakonometrics.free.fr/US.txt", 
+ header=TRUE,sep=";")

(je renvoie sur le précédant billet pour un descriptif précis). Idéalement, ça serait parfait si on pouvait visualiser sur une carte les variables. Pour cela, il faut rajouter une colonne à notre base, avec le nom complet des états,

> abreviation=read.table( 
+ "http://freakonometrics.blog.free.fr/public/data/etatus.csv", 
> header=TRUE,sep=",") 
> US$USPS=rownames(US) 
> US=merge(US,abreviation) 
> US$nom=tolower(US$NOM)

Cette fois, on va pouvoir faire de la cartographie, les noms de nos états étant (presque) les mêmes que ceux des cartes de R,

> library(maps) 
> VL0=strsplit(map("state")$names,":") 
> VL=VL0[[1]] 
> for(i in 2:length(VL0)){VL=c(VL,VL0[[i]][1])} 
> ETAT=match(VL,US$nom)

Cette fois-ci, on a toutes les informations pour faire une carte, avec une couleur fonction de la variable d’intérêt (espérance de vie à la naissance, taux d’homicides, taux illettrisme, etc).

> library(RColorBrewer) 
> carte=function(V=US$Murder,titre= 
+ "Taux d'homicides aux Etats-Unis"){ 
+ variable=as.numeric(as.character(cut(V, 
+ quantile(V,seq(0,1,by=1/6)),labels=1:6))) 
+ niveau=variable[ETAT] 
+ couleur=rev(brewer.pal(6,"RdBu")) 
+ noml=levels(cut(V,quantile(V,seq(0,1,by=1/6)))) 
+ map("state", fill = TRUE, col=couleur[niveau]); 
+ legend(-78,34,legend=noml,fill=couleur, + cex=1,bty="n"); 
+ title(titre)}

Commençons par l’analyse du nombre de jours durant lesquels la température passe en dessous de 0°C, par an, en moyenne, dans la capitale (ou la plus grande ville) de l’état, afin de tester notre fonction,

> carte(US$Frost, + titre="Nombre de jours de gel par an")

Pour le taux d’homicide (qui est la variable par défaut) on a

> carte()

Ça sera notre variable d’intérêt lors de la modélisation de mardi. Enfin, on peut lancer une régression, et représenter spatialement les résidus,

> reg=lm(Murder~.-NOM-USPS-nom,data=US) 
> regs=step(reg) 
> carte(residuals(regs), 
+ titre="Résidus de la régression")

Nous voila équipés pour commencer l’économétrie spatiale…

Pour aller maintenant un peu plus loin dans la modélisation, je vais rajouter une variable qualitative, par exemple l’appartenance politique du gouverneur en 1977 (les données datent de cette époque). Les données ont été extraites de wikipedia, suite aux élections de 1974, 1975, 1976 et 1977,

> GV=read.table( 
+ "http://freakonometrics.blog.free.fr/public/data/governor.csv", 
+ header=TRUE,sep=";") 
> etat=strsplit(as.character(GV$State),"-") 
> listeetat=rep(NA,nrow(GV)) 
> for(i in 1:nrow(GV)){ 
+ listeetat[i]=etat[[i]][1] 
+ } 
> indice=which(is.na(listeetat)==FALSE) 
> basegv=data.frame(state=tolower(listeetat[indice]), 
+ party=GV$Party[indice])

On a la visualisation suivante

> library(maps) 
> library(RColorBrewer) 
> couleur=rev(brewer.pal(6, "RdBu")) 
> Z=rep(6,length(basegv$party)) 
> Z[basegv$party=="Democratic"]=1 
> VL0=strsplit(map("state")$names,":")
> VL=VL0[[1]] 
> for(i in 2:length(VL0)){VL=c(VL,VL0[[i]][1])} 
> ETAT=match(VL,basegv$state)
> niveau=Z[ETAT] 
> map("state", fill = TRUE, col=couleur[niveau])

For those who think more variates should be added to the dataset, some can be found e.g. on http://www.statemaster.com/, like the total executions since 1930, or the date the state joint the U.S.A.

Régression, variables explicatives et géométrie

La régression (comme tout calcul d’espérance conditionnelle) est un problème de projections

Un théorème intéressant est le théorème dit de Frisch-Waugh, permettant de comprendre la différence fondamentale entre un modèle de régression multiple, et les modèles de régression simple. On reviendra sur ce point en cours en évoquant rapidement le paradoxe de Simpson. La formulation est la suivante: on veut régresser https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW15.gif sur https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW16.gif et https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW17.gif, deux ensembles (a priori disjoints) de variables explicatives,


Les équations normales (associées au problème de minimisation de la somme des carrés des erreurs) sont


de telle sorte qu’à l’optimum on peut relier les estimateurs des deux jeux de paramètres par


La première partie correspond à la régression de https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW15.gif sur https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW13.gif, mais il reste un second terme dès lors que https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW16.gif et https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW17.gif ne sont pas orthogonales. Notons https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW04.gifla matrice de projection (orthogonale) sur https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW16.gif


Alors on peut écrire


En posant

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW07.gif et https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW08.gif,

on retrouve un modèle linéaire classique (à condition de travailler sur la projection des variables sur le sous-espace engendré par https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/05/FW16.gif , i.e. en transformant les variables),


Pour aller plus loin sur la géométrie des moindres carrées, et sur le théorème de Frisch-Waugh, je peux renvoyer à des notes de cours, et à quelques transparents

Le graphique ci-dessous correspond au cas où les variables explicatives sont orthogonales, et dans ce cas, la régression multiple est équivalente à deux régressions simples


Le graphique ci-dessous au cas non orthogonal,


Le code qui permet de vérifier ces histoires de projections successives est relativement simple. Tout d’abord on importe les données, et on regarde le modèle global,

> chicago=read.table(
+ "http://freakonometrics.free.fr/chicago.txt",
+ header=TRUE,sep=";")
> Y=chicago$Fire
> X1=chicago$X_1
> X2=chicago$X_2
> X3=chicago$X_3
> base=data.frame(Y,X1,X2,X3)
> tail(base)
Y    X1 X2     X3
42  4.8 0.152 19 13.323
43 10.4 0.408 25 12.960
44 15.6 0.578 28 11.260
45  7.0 0.114  3 10.080
46  7.1 0.492 23 11.428
47  4.9 0.466 27 13.731
> regression=lm(Y~X1+X2+X3)
> summary(regression)

lm(formula = Y ~ X1 + X2 + X3)

Min     1Q Median     3Q    Max
-9.737 -4.565 -1.479  3.751 16.079

Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.07525    6.19447   3.564 0.000910 ***
X1          -0.62764    5.28130  -0.119 0.905953
X2           0.22378    0.06161   3.632 0.000744 ***
X3          -1.55059    0.38195  -4.060 0.000204 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.527 on 43 degrees of freedom
Multiple R-squared: 0.4417,	Adjusted R-squared: 0.4027
F-statistic: 11.34 on 3 and 43 DF,  p-value: 1.314e-05

> n=length(Y)
> X=matrix(c(rep(1,n),X1,X2,X3),n,4)
> X[1:5,]
[,1]  [,2] [,3]   [,4]
[1,]    1 0.604   29 11.744
[2,]    1 0.765   44  9.323
[3,]    1 0.735   36  9.948
[4,]    1 0.669   37 10.656
[5,]    1 0.814   53  9.730

Ensuite, on va projeter seulement sur les deux premières variables (et la constante)

> FWX1=X[,1:3]
> regression12=lm(Y~X1+X2)
> solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%Y
[1,]  0.08069764
[2,] 11.56913900
[3,]  0.15108490
> summary(regression12)$coefficients
Estimate Std. Error    t value   Pr(>|t|)
(Intercept)  0.08069764 3.49182154 0.02311047 0.98166664
X1          11.56913900 5.05014782 2.29085156 0.02681947
X2           0.15108490 0.06854408 2.20420043 0.03278652
> FWX2=X[,4]
> H1=FWX1%*%solve(t(FWX1)%*%FWX1)%*%t(FWX1)
> M1=diag(rep(1,n))-H1
> FWX2s=M1%*%FWX2
> FWYs =M1%*%Y
> (beta2=solve(t(FWX2s)%*%FWX2s)%*%t(FWX2s)%*%FWYs)
[1,] -1.550594
> summary(regression)$coefficients[4]
[1] -1.550594
> (beta1=solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%Y-
+        solve(t(FWX1)%*%FWX1)%*%t(FWX1)%*%FWX2%*%beta2)
[1,] 22.0752495
[2,] -0.6276442
[3,]  0.2237765
> summary(regression)$coefficients[1:3]
[1] 22.0752495 -0.6276442  0.2237765

On retrouve bien l’estimateur du dernier paramètre à l’aide du théorème de Frisch-Waugh, puis en utilisant les équations normales, on en déduit les trois premiers.

Tail index estimation

These data were collected at Copenhagen Reinsurance and comprise 2167 fire losses over the period 1980 to 1990, They have been adjusted for inflation to reflect 1985 values and are expressed in millions of Danish Kron. Note that it is possible to work with the same data as above but the total claim has been divided into a building loss, a loss of contents and a loss of profits.

> base1=read.table(
+ "http://freakonometrics.free.fr/danish-univariate.txt",
+ header=TRUE)
> base2=read.table(
+ "http://freakonometrics.free.fr/danish-multivariate.txt",
+ header=TRUE)

Consider here the first dataset (we deal – so far – with univariate extremes),

> X=base1$Loss.in.DKM
> D=as.Date(as.character(base1$Date),"%m/%d/%Y")
> plot(D,X,type="h")

The graph is the following,

A natural idea is then to plot



> Xs=sort(X)
> logXs=rev(log(Xs))
> n=length(X)
> plot(log(Xs),log((n:1)/(n+1)))

Points are on a straight line here. The slope can be obtained using a linear regression,

> B=data.frame(X=log(Xs),Y=log((n:1)/(n+1)))
> reg=lm(Y~X,data=B)
> summary(reg)

lm(formula = Y ~ X, data = B)

Min       1Q   Median       3Q      Max
-0.59999 -0.00777  0.00878  0.02461  0.20309

Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.089442   0.001572   56.88   <2e-16 ***
X           -1.382181   0.001477 -935.55   <2e-16 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.04928 on 2165 degrees of freedom
Multiple R-squared: 0.9975,	Adjusted R-squared: 0.9975
F-statistic: 8.753e+05 on 1 and 2165 DF,  p-value: < 2.2e-16

> reg=lm(Y~X,data=B[(n-500):n,])
> summary(reg)

lm(formula = Y ~ X, data = B[(n - 500):n, ])

Min       1Q   Median       3Q      Max
-0.48502 -0.02148 -0.00900  0.01626  0.35798

Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.186188   0.010033   18.56   <2e-16 ***
X           -1.432767   0.005105 -280.68   <2e-16 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.07751 on 499 degrees of freedom
Multiple R-squared: 0.9937,	Adjusted R-squared: 0.9937
F-statistic: 7.878e+04 on 1 and 499 DF,  p-value: < 2.2e-16

> reg=lm(Y~X,data=B[(n-100):n,])
> summary(reg)

lm(formula = Y ~ X, data = B[(n - 100):n, ])

Min       1Q   Median       3Q      Max
-0.33396 -0.03743  0.02279  0.04754  0.62946

Estimate Std. Error t value Pr(>|t|)
(Intercept)  0.67377    0.06777   9.942   <2e-16 ***
X           -1.58536    0.02240 -70.772   <2e-16 ***
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.1299 on 99 degrees of freedom
Multiple R-squared: 0.9806,	Adjusted R-squared: 0.9804
F-statistic:  5009 on 1 and 99 DF,  p-value: < 2.2e-16

The slope here is somehow related to the tail index of the distribution. Consider some heavy tailed distribution, i.e. https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill03.gif, so that https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill27.gif, where https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill28.gif is some slowly varying function. Equivalently, the exists a slowly varying function https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill29.gif such that https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill30.gif. Then


i.e. since a natural estimator for https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill35.gif is the order statistic https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill36.gif, the slope of the straight line is the opposite of tail index https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill98.gif. The estimator of the slope is (considering only the https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill99.gif largest observations)


Hill‘s estimator is based on the assumption that the denominator above is almost 1 (which means that  https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill15.gif, as https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill16.gif), i.e.


Note that, if https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill14.gif, but not two fast, i.e. https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill15.gif as https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill16.gif, then https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill12.gif (one can even get https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill11.gif  with stronger convergence assumptions). Further


Based on that (asymptotic) distribution, it is possible to get a (asymptotic) confidence interval for https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill98.gif

> xi=1/(1:n)*cumsum(logXs)-logXs
> xise=1.96/sqrt(1:n)*xi
> plot(1:n,xi,type="l",ylim=range(c(xi+xise,xi-xise)),
+ xlab="",ylab="",)
> polygon(c(1:n,n:1),c(xi+xise,rev(xi-xise)),
+ border=NA,col="lightblue")
> lines(1:n,xi+xise,col="red",lwd=1.5)
> lines(1:n,xi-xise,col="red",lwd=1.5)
> lines(1:n,xi,lwd=1.5)
> abline(h=0,col="grey")

It is also possible to work with https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill06.gif, then https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill05.gif. And similarly https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill13.gif as https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill14.gif (and again https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill10.gif with additional assumptions on the rate of convergence), and


(obtained using the delta-method). Again, we can use that result to derive (asymptotic) confidence intervals

> alpha=1/xi
> alphase=1.96/sqrt(1:n)/xi
> YL=c(0,3)
> plot(1:n,alpha,type="l",ylim=YL,xlab="",ylab="",)
> polygon(c(1:n,n:1),c(alpha+alphase,rev(alpha-alphase)),
+ border=NA,col="lightblue")
> lines(1:n,alpha+alphase,col="red",lwd=1.5)
> lines(1:n,alpha-alphase,col="red",lwd=1.5)
> lines(1:n,alpha,lwd=1.5)
> abline(h=0,col="grey")

The Deckers-Einmahl-de Haan estimator is


where for


Then (given again conditions on the speed of convergence i.e. https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill14.gif, with https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill15.gif as https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill16.gif),


Finally, Pickands‘ estimator


it is possible to prove that, as https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/hill14.gif,


Here the code is

> Xs=rev(sort(X))
> xi=1/log(2)*log( (Xs[seq(1,length=trunc(n/4),by=1)]-
+ Xs[seq(2,length=trunc(n/4),by=2)])/
+ (Xs[seq(2,length=trunc(n/4),by=2)]-Xs[seq(4,
+ length=trunc(n/4),by=4)]) )
> xise=1.96/sqrt(seq(1,length=trunc(n/4),by=1))*
+sqrt( xi^2*(2^(xi+1)+1)/((2*(2^xi-1)*log(2))^2))
> plot(seq(1,length=trunc(n/4),by=1),xi,type="l",
+ ylim=c(0,3),xlab="",ylab="",)
> polygon(c(seq(1,length=trunc(n/4),by=1),rev(seq(1,
+ length=trunc(n/4),by=1))),c(xi+xise,rev(xi-xise)),
+ border=NA,col="lightblue")
> lines(seq(1,length=trunc(n/4),by=1),
+ xi+xise,col="red",lwd=1.5)
> lines(seq(1,length=trunc(n/4),by=1),
+ xi-xise,col="red",lwd=1.5)
> lines(seq(1,length=trunc(n/4),by=1),xi,lwd=1.5)
> abline(h=0,col="grey")

It is also possible to use maximum likelihood techniques to fit a GPD distribution over a high threshold.

> library(evd)
> library(evir)
> gpd(X,5)
[1] 2167

[1] 5

[1] 0.8827873

[1] 254

[1] "ml"

xi      beta
0.6320499 3.8074817

xi      beta
0.1117143 0.4637270

[,1]        [,2]
[1,]  0.01248007 -0.03203283
[2,] -0.03203283  0.21504269

[1] "observed"

[1] 0

[1] 754.1115

[1] "gpd"

or equivalently (or almost)

> gpd.fit(X,5)
[1] 5

[1] 254

[1] 0

[1] 754.1115

[1] 3.8078632 0.6315749

[1] 0.1172127

[1] 0.4636270 0.1116136

The interest of the latest function is that it is possible to visualize the profile likelihood of the tail index,

> gpd.profxi(gpd.fit(X,5),xlow=0,xup=3)


> gpd.profxi(gpd.fit(X,20),xlow=0,xup=3)

Hence, it is possible to plot the maximum likelihood estimator of the tail index, as a function of the threshold (including a confidence interval),

> GPDE=Vectorize(function(u){gpd(X,u)$par.ests[1]})
> GPDS=Vectorize(function(u){
+ gpd(X,u)$par.ses[1]})
> u=c(seq(2,10,by=.5),seq(11,25))
> XI=GPDE(u)
> plot(u,XI,ylim=c(0,2))
> segments(u,XI-1.96*XIS,u,XI+
+ 1.96*XIS,lwd=2,col="red")

Finally, it is possible to use block-maxima techniques.

> gev.fit(X)
[1] 0

[1] 3392.418

[1] 1.4833484 0.5930190 0.9168128

[1] 0.01507776 0.01866719 0.03035380

The estimator of the tail index is here the last coefficient, on the right.
Since it is rather difficult to install a package in class rooms, here is the source of rcodes used here (to fit a GPD for exceedances)

> source("http://freakonometrics.blog.free.fr/public/code/gpd.R")

Next time, we will discuss how to use those estimators.

Quel lancer de dé faire à Serpents et Échelles ?

Hier soir, j’évoquais l’utilisation des chaînes de Markov au jeu Serpents et Échelles. Comme me le faisait remarquer Jean-Philippe, de manière assez étrange, les enfants sont toujours beaucoup plus content après avoir fait un 6 qu’après avoir fait un 1. Mais est-ce réellement la valeur optimale que le dé doit prendre ? Ça dépend bien sûr de la position. Par exemple, au premier lancer, on peut se demander ce que deviendrait le nombre de tours moyen nécessaires pour finir le jeu, conditionnellement aux 6 lancers possibles. On peut calculer, toujours conditionnellement aux 6 valeurs possibles de dés (et des positions où on va se retrouver), l’espérance du nombre de tours pour finir la partie,

INITIAL = as.numeric(which(M[h0+1,]>0))-1
for(k in 1:length(INITIAL)){
initial=rep(0,n+1); initial[INITIAL[k]]=1
for(h in 1:length(game)){

(à partir du code mis en ligne hier), e.g. pour le premier lancer,

> esperance(0)
1        2        3        5        6       14
32.16499 31.99947 31.82348 31.47954 31.36441 29.83020

où la valeur indiquée est la position où l’on pourrait se retrouver (la dernière correspond à un 4). On note que le meilleur premier lancer possible est celui qui nous amène sur la première échelle, i.e. la 4ème case,

Si on regarde case par case, on a les valeurs des “meilleurs” lancers de dés suivant (sans forcément l’unicité) en fonction de sa position sur le plateau

(avec en rouge les serpents et en bleu les échelles). Amusant non ?

Wald, score et rapport de vraisemblance

Vendredi, nous devrions aborder un peu la problématique des tests asymptotiques. Avant, rappelons que l’estimateur du maximum de vraisemblance de  vérifie la propriété asymptotique suivante.

Or on a une propriété intéressante sur les convergences de vecteurs Gaussiens (pour rester assez général). Soit  une suite de vecteurs aléatoires telle que  où . Alors

où  (pour un paramètre univarié, ça sera alors ). Mais ce résultat n’est intéressant que si  est connue. En fait, si on a la suite des variances vérifie , alors

Traduit sur notre estimateur du maximum de vraisemblance, cela signifie que

ou encore

Trois tests peut alors être mis en œuvre (et tant qu’à faire, on peut essayer de l’illustrer sur un exemple: encore et toujours du pile/face).

> set.seed(1)
> n=20
> X=sample(0:1,size=n,replace=TRUE)
> p=seq(0,1,by=.01)
> logL=function(p){sum(log(dbinom(X,size=1,prob=p)))}
> LL=sapply(p,logL)
> plot(p,LL,type="l",col="red",lwd=2)
> p0=.5
> points(p0,logL(p0),pch=3,cex=1.5,lwd=2)
> abline(v=p0,lty=2)

On a alors 20 tirages de pile ou face, on a obtenu 11 piles, et on se demande si la pièce est “juste“.Mais avant toute chose, commençons par calculer l’estimateur du maximum de vraisemblance, ainsi que le score et l’information de Fisher. Dans ce modèle de variables de Bernoulli, on connaît des formes explicites de ces quantités. Mais ici, on va utiliser la méthode la plus générale qui soit, i.e. on va maximiser la log-vraisemblance, puis on va calculer le score et l’information de Fisher à la valeur de  que l’on souhaite tester.

> neglogL=function(p){-sum(log(dbinom(X,size=1,prob=p)))}
> pml=optim(fn=neglogL,par=p0,method="BFGS")$par
Warning messages:
1: In dbinom(x, size, prob, log) : NaNs produced
2: In dbinom(x, size, prob, log) : NaNs produced
> nx=sum(X==1)
> f = expression(nx*log(p)+(n-nx)*log(1-p))
> Df = D(f, "p")
> Df2 = D(Df, "p")
> p=p0
> score=eval(Df)
> (IF=-eval(Df2))
[1] 80
> 1/(p0*(1-p0)/n)
[1] 80

On note que l’information de Fisher calculée par double différenciation de la log-vraisemblance est identique à la formule analytique. Pour rappels, on a

Trois tests (équivalent asymptotiquement comme on le verra en cours) sont alors possibles.

Tout d’abord le test de Wald propose de travailler sur la différence entre l’estimateur du maximum de vraisemblance, et la valeur que l’on cherche à tester (comme le montre le dessin ci-dessous). Cette différence doit être “petite” si  est la vraie valeur.On peut alors utiliser la statistique suivante

Asymptotiquement, en utilisant le résultat évoqué au début, cette statistique tend, sous l’hypothèse que  est la vraie valeur, vers une loi du chi-deux.

> alpha=0.05
> (T=(pml-p0)^2*IF)
[1] 0.1999970
> T<qchisq(1-alpha,df=1)
[1] TRUE

Ici, on accepte l’hypothèse que la pièce n’est pas pipée.

Ensuite le test du rapport de vraisemblance propose de travailler sur les valeurs de la log-vraisemblance. Là encore, cette différence doit être “petite” si  est la vraie valeur. Graphiquement, on mesure une distance normalisée

La statistique est celle d’un test bilatéral,

On pose alors

(le 2 sera détaillé un peu en cours, c’est mis là de manière à avoir trois tests équivalents). Asymptotiquement, on peut montrer que cette statistique tend, sous l’hypothèse que  est la vraie valeur, vers une loi du chi-deux.

> (T=2*(logL(pml)-logL(p0)))
[1] 0.2003347
> T<qchisq(1-alpha,df=1)
[1] TRUE

Là encore, on accepte l’hypothèse que la pièce n’est pas pipée.

Enfin, le test du score propose de travailler sur la pente de la log-vraisemblance en  . Cette pente doit être “petite” si  est la vraie valeur.

La pente correspond au score,

La statistique de test est alors ici

Et là encore, asymptotiquement, on peut montrer que cette statistique tend, sous l’hypothèse que  est la vraie valeur, vers une loi du chi-deux.

> (T=slope^2/IF)
[1] 0.2
> T<qchisq(1-alpha,df=1)
[1] TRUE

Et le test accepte ici encore l’hypothèse que la pièce n’est pas pipée.

Extracting information from a keyboard…

Yesterday, Baptiste published a post on “ethno-photography” (here). As he mentioned it, in Paris 8, they experience a real absence of serious cleaning of office equipment. He then shows the keyboard of the only computer they can use in the sociology department (for forty researchers),

Apart from the fact that everyone in France should be ashamed to see how much is spent in universities (which is the first information we have from that picture), we should also be able to guess in which langage people work in this department.
I considered three books (two in French, one in English) and I would like to see the frequency of each letter,

  • Mauss, manuel d’éthnographie (here), 1926
  • Durkheim, Livre II: Les croyances élémentaires in Les formes élémentaires de la vie religieuse (here), 1912
  • Ferri, Criminal Sociology (here), 1896

Those three books are in rich text format, I just changed it to get text files… Then, it is easy to count appearance of letters. E.g. for Mauss,

> library(corpora)
> textfile=scan("MAUSS-manuel.txt",
+ what="char", sep="\n")
Read 1550 items
> textfile<-tolower(textfile)
> M=NA
> for(i in 1:length(textfile)){
+ line=textfile[i]
+ M=c(M,strsplit(as.character(line),"")[[1]])
+ }
> T=table(M)
> T
    '     -           \t     !     "     %     &     (     )     ,     .     / 
 5308  1049 86589    44     3     3     2     2   370   391  6609  4909    12 
    :     ;     ?     @     ]     _     ~     ’     =     «     »     ¬     ° 
  819  1178   113     1     1     4     1    39     1   108   107   823     3 
    …     0     1     2     3     4     5     6     7     8     9     a     à 
    1    69   213    83    73    34    48    33    28    64   151 30559  1651 
    â     ä     b     c     ç     d     e     é     è     ê     ë     f     g 
  224     3  3562 14678   110 17713 63955 10354  1798  1000     5  4555  4911 
    h     i     î     ï     j     k     l     m     n     ñ     o     º     ô 
 4359 30851   226    47  1147   247 24792 12844 32525     6 25562     2   151 
    ö     œ     p     q     r     s     t     u     ù     û     ü     v     w 
   12    52 12696  4667 28237 37630 32945 25001   211    40     9  4787   164 
    x     y     z 
 1996  1222   343

Then, we can summarize in to see proportion of standard 26 letters, and we have, for Mauss,

and for Durkeim,

If we compare the two, we have almost the same proportions,

If we look at our book in English now, we have

i.e., if we compare with Mauss for instance

So we have much more E in French than in English, but still, people writing in English use a lot the E. So looking at the E should not give us any clue…. But we can see that in English, the H is as common as the L, or the C. Not in French, where L is much more frequent than the H. But on the picture, the C is more clear than the H. We can also look at the U, which is common in French, not in English… Here, on the keyboard, it is perfectly clear… so I guess people use it frequently.
So I would say that they write more in French than the write in English, on that computer.
Actually, the same idea has been used a long time ago on calculators to see that Benford’s law works: some numbers are really used (as well as the legend pretends that some pages in logarithm books were never used….), see here orthere. So Baptiste, if one day the keyboard is cleaned up, please send me another picture after a few weeks to see if things have changed….
An for those who cannot imagine how it is to work in some universities in France, just look at his blog (here). Pictures are unbelievable….Good luck Baptiste….

Qui a eu l’idée d’élargir les trottoirs en France ?

J’étais tombé l’autre jour sur un joli dessin expliquant pourquoi les policiers et les syndicats n’estimaient jamais le même nombre de manifestants, avec une explication scientifique, et vaguement rationnelle, ici.

C’est  effectivement joli… et on retrouve un rapport de 3 dont pas mal de monde semble parler, sauf que si c’était vrai, le rapport entre les estimations de la police et des syndicats devrait être stable d’une manifestation à l’autre, non ?
Si on compare le 19 mars dernier (ici) et le 16 octobre (), on arrive aux chiffres suivants,

Les traits pointillés sont la droite de régression. La courbe en trait plein et la régression (passant par l’origine car on cherche un rapport de proportionnalité) en prenant en compte un effet taille, i.e. en introduisant des poids différents dans ma régression. J’ai pris la taille des villes, autrement dit, plus la ville est grosse, plus je donne un poids important au rapport… L’idée étant d’intégrer un effet médiatique, même si la pente n’est plus le ratio moyen…  Avec cette dernière méthode, j’obtiens des pentes de 4,06 pour mars, et 6,12 pour octobre… ce qui est à un peu différents des rapports nationaux mentionnés sur le site du monde, i.e. (3 millions selon les syndicats et 1,2 selon la police, en mars, contre  2,5 millions contre 825 milliers l’autre jour). Notons qu’on peut vérifier rapidement que ce n’est pas un problème statistique, les valeurs sont sensiblement différentes, si l’on compare les distributions des estimateurs..

Bref, quelque chose semble avoir changé entre ces deux dates…. la taille des trottoirs ?

Margin of error, and comparing proportions in the same sample

Irecently tried to answer a simple question, asked by @adelaigue. Actually, I thought that the answer would be obvious… but it is a little bit more compexe than what I thought. In a recent survey about elections in Brazil, it was mentionned in a French newspapper that “Mme Rousseff, 62 ans, de 46,8% des intentions de vote et José Serra, 68 ans, de 42,7%” (i.e. proportions obtained from the survey). It is also mentioned that “la marge d’erreur du sondage est de 2,2% ” i.e. the margin of error is 2.2%, which means (for the journalist) that there is a “grande probabilité que les 2 candidats soient à égalité” (there is a “large probability” to have equal proportions).
Usually, in sampling theory, we look at the margin of error of a single proportion. The idea is that the variance of https://latex.codecogs.com/gif.latex?%20\widehat{p}, obtained from a sample of size https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial15.png is


thus, the standard error is


The standard 95% confidence interval, derived from a Gaussian approximation of the binomial distribution is


The largest value is obtained when p is 1/2, and then we have a worst case confidence interval (an upper bound) which is


So with a margin of error https://perso.univ-rennes1.fr/arthur.charpentier/latex/m205.png means that https://perso.univ-rennes1.fr/arthur.charpentier/latex/m206.png. Hence, with a 5% margin of error, it means that n=400. While 2.2% means that n=2000:
> 1/.022^2
[1] 2066.116
Classically, we compare proportions between two samples: surveys at two different dates, surveys in different regions, surveys paid by two different newpapers, etc. But here, we wish to compare proportions within the same sample. This has been consider in an “old” paper published in 1993 in the American Statistician,

It contains nice figures to illustrate the difference between the standard approach,

and the one we would like to study here.

This point is mentioned in the book by Kish, survey sampling (thanks Benoit for the reference),

Let https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial05.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial06.png denote empirical frequencies we have obtained from the sample, based on https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial15.png observations. Then since




we have


Thus, a natural margin of error on the difference between the two proportion is here


which is here 4 points
> n=2000
> p1=46.8/100
> p2=42.7/100
> 1.96*sqrt((p1+p2)-(p1-p2)^2)/sqrt(n)
[1] 0.04142327
Which is exactly the difference we have here ! Hence, the probability of reaching such a value is quite small (2%)
> s=sqrt(p1*(1-p1)/n+p2*(1-p2)/n+2*p1*p2/n)
> (p1-p2)/s
[1] 1.939972
> 1-pnorm(p1-p2,mean=0,sd=sqrt((p1+p2)-(p1-p2)^2)/sqrt(n))
[1] 0.02619152

Actually, we can compare the three margin of errors we have so far,

  • the upper bound
  • the “average one”


  • the more accurate one we just obtained,

where https://perso.univ-rennes1.fr/arthur.charpentier/latex/m214.png.
> p=seq(0,.5,by=.01)
> ic1=rep(1.96/sqrt(4*n),length(p))
> ic2=1.96*sqrt(p*(1-p))/sqrt(n)
> delta=.01
> ic31=1.96*sqrt(2*p-delta^2)/sqrt(n)
> delta=.2
> ic32=1.96*sqrt(2*p-delta^2)/sqrt(n)
> plot(p,ic32,type=”l”,col=”blue”)
> lines(p,ic31,col=”red”)
> lines(p,ic2)
> lines(p,ic1,lty=2)
So on the graph below, the dotted line is the standard upper bound, the plain line in black being a more accurate one when the probability is https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial19.png (the x-axis). The red line is the true margin of error with a large difference between candidates (20 points) and the blue line with a small difference (1 point).

Remark: an alternative is to consider a chi-square test, comparering two multinomial distributions, with probabilities https://perso.univ-rennes1.fr/arthur.charpentier/latex/m215.png and https://perso.univ-rennes1.fr/arthur.charpentier/latex/m216.png where https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial19.png is the average proportion, i.e. 44.75%. Then


i.e.  https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial22.png=3.71
> p=(p1+p2)/2
> (x2=n*((p1-p)^2/p+(p2-p)^2/p))
[1] 3.756425
> 1-pchisq(x2,df=1)
[1] 0.05260495
Under the null hypothesis, https://perso.univ-rennes1.fr/arthur.charpentier/latex/multinomial22.png should have a chi-square distribution, with one degree of freedom (since the average is fixed here). Here the probability to reach that level is around 5% (which can be compared with the 2% we add before).

So finally, I would think that here, stating that there is a “large probability” is not correct…

Too large datasets for regression ? What about subsampling….

recently, a classmate working in an insurance company told me he had too large datasets to run simple regressions (GLM, which involves optimization issues), and that they were thinking of a reward for the one who will write the best R-code (at least the fastest). My first idea was to use subsampling techniques, saying that 10 regressions on 100,000 observations can take less time than a regression on 1,000,000 observations. And perhaps provide also better results…

  • Time to run a regression, as a function of the number of observations

Here, I generate a dataset as follows


and we fit


where https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp03.png is a spline function (just to make it as general as possible, since in insurance ratemaking, we include continuous variates that do not influence claims frequency linearly in the score). Yes, there might be also useless variables, including one of them which is strongly correlated with one that has an impact in the regression. The code to generate the dataset is simply

> n=10000
> X1=rexp(n)
> X2=sample(c("A","B","C"),size=n,replace=TRUE)
> X3=runif(n)
> Z=rmnorm(n,c(0,0),matrix(c(1,0.8,.8,1),2,2))
> X4=Z[,1]
> X5=Z[,2]
> X6=X1^2
> E=runif(n)
> lambda=.2*X5-4*dbeta(X3,2,5)+X1+
> Y=rpois(n,exp(lambda))
> base=data.frame(Y,X1,X2,X3,X4,X5,X6,E)

We would like the study the time it takes to run a regression, as a function of the size (i.e. the number of lines https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp04.png) of the dataset.

> system.time( glm(Y~bs(X1)+X2+X3+X4+
+ X5+X6+offset(log(E)),family=poisson,
+ data=base) )
utilisateur     système      écoulé
0.25        0.00        0.25

Here, the time I look at is the last one. But so far, it was rather simple, but it is not the best model I can get. Let us use a stepwise (backward) variable selection,

> system.time( step(glm(Y~bs(X1)+X2+X3+
+ X4+X5+X6+offset(log(E)),family=poisson,
+ data=base)) )
Start:  AIC=2882.1
Y ~ bs(X1) + X2 + X3 + X4 + X5 + X6 + offset(log(E))
Step:  AIC=2882.1
Y ~ bs(X1) + X2 + X3 + X4 + X5 + offset(log(E))
Df Deviance    AIC
<none>        2236.0 2882.1
- X5      1   2240.1 2884.2
- X4      1   2244.1 2888.2
- X3      1   4783.2 5427.3
- X2      2   5311.4 5953.5
- bs(X1)  3   6273.7 6913.8
utilisateur     système      écoulé
1.82        0.03        1.86

Finally, from the first regression, we have points in black (based on 200 simulated datasets), and with a stepwise procedure, we have the points in red.

i.e. it might look linear (proportional), but if it was linear, then on a log-log scale, we should have also straigh lines, with slope 1,

Actually, it looks like a convex function.

The interpretation of that convexity might lead to misinterpretation. On the graph below on the left, on a dataset two times bigger than the previous one (black point) will be less than two times longer to run, while on the right, it will be more than two timess longer,

Convexity can simply be interpreted as “too large datasets take time, and too small too…”. Which is a first step: it should be interesting, in some cases, to run several regressions on smaller datasets….

  • Running 100 regressions on 100 lines, or running 1 regression on 10,000 lines ?

Here, we have datasets with https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp04.png=200,000 lines. The questions is how long will it take if we subdived into https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp05.png subsamples (of equal size), and run https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp05.png regressions ?

> nk=trunc(n/k)rep(1:k,each=nk); nt=nk*k
> base=data.frame(Y[1:nt],X1[1:nt],
+ X2[1:nt],X3[1:nt],X4[1:nt],X5[1:nt],
+ X6[1:nt],E[1:nt],classe)
> system.time( for(j in 1:k){
+  glm(Y~bs(X1)+X2+X3+X4+X5+
+ X6+offset(log(E)),family=poisson
+ ,data=base,subset=classe==j) })
utilisateur     système      écoulé
1.31        0.00        1.31
> system.time( for(j in 1:k){
+      step(glm(Y~bs(X1)+X2+X3+
+ X4+X5+X6+offset(log(E)),family=
+ poisson,data=base,subset=classe==j)) })
Start:  AIC=183.97
Y ~ bs(X1) + X2 + X3 + X4 + X5 + X6 + offset(log(E))


  Df Deviance    AIC
<none>        117.15 213.04
- X2      2   250.15 342.04
- X3      1   251.00 344.89
- X4      1   420.63 514.53
- bs(X1)  3   626.84 716.74
utilisateur     système      écoulé
11.97        0.03       12.31

On the graph below, we have the time (y-axis, here on a log scale) it took to run https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp05.png regression on samples of size https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp06.png, as function of https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp05.png (x-axis), including the time it took to run the regression on a dataset of size https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp04.png which is the concentration of dots on the left (i.e. https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp05.png=1), both on the 6 regressors – in black – and with a strepwise procedure – in red. One has to keep in mind that I did not remove the printing option in the stepwise procedure, so it might be difficult to compare the two clouds (black vs. red). Nevertheless, we clearly see that if we run https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp05.png regression on samples of size https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp06.png, when https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp05.png is not too large, i.e. less than 10 or 15, it is not longer than the regression on https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp04.png=200,000 lines.

So here we see that running 100 regressions on 2,000 lines is longer than running 1 regression on 200,000 lines… But maybe we are not comparing things that are actually comparable: what if it takes a bit longer, but we strongely improve the quality of our estimators ?

  • What about the quality of the output ?

Here, we consider only one dataset, with https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp04.png=100,000 lines (just to make it run a bit faster). And https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp05.png=20 subsets. Recall that the generated dataset is from


and we fit


Here, we plot here https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp07.png and a confidence interval, defined as


The lightblue segment is the initial estimator, while the blue one is obtained from the stepwise procedure. The grey area represent the estimation on the overall sample, while the https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp05.png segments on the right are the https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp05.png estimators (each on samples of size https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp06.png).

We can see that we have much more volatility on those https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp05.png estimators, but the average (horizontal doted lines) are not so bad… The true value (i.e. the one used to generate the dataset is the dotter black horizontal line).
And if we repeat that on 1,000 simulated dataset, we obtaind the following distribution for https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp07.png (blue line), so we have an unbiased estimator of our parameter (the verticular line being here the true value), here including a stepwise procedure,

But if we add the the red curve is the average of the https://f.hypotheses.org/wp-content/blogs.dir/253/files/2016/11/largesamp09.png the previous one being now the clear blue line in the back, we see that taking average of estimators on subsamples is not bad at all, on the contrary,

and for those who think that the stepwise procedure is a mistake, here is what we get without it,

So what we can see is that running 20 regressions can take (a little) more time (from what we’ve seen earlier) than running only one on the whole dataset…. but it provides better estimates. So the tradeoff is not that simple, and maybe running several regressions on huge datasets can be a proper alternative.

Tails of copulas, une lecture graphique

Suite à une formation que je faisais en fin de semaine à Brest (les slides sont ici et ), je voulais revenir sur les histoires de tails of copulas, pour reprendre le titre de l’article (ici) de Gary Venter (et qui correspond à des choses que j’avais pu présenter il y a quelques années à Berlin, les slides étant en ligne ici).

  • Quantifier la dépendance de queue

L’idée est de noter qu’il est noter qu’il existe deux manières de quantifier la dépendance de queue. La première est liée à l’approche de Joe (1990, ici, ou 1997 pour le livre), qui a introduit un (strong) tail dependence index. Par exemple pour la queue inférieure,




La seconde est liée à une idée que l’on retrouve dans les travaux de Janet Heffernan, Stuart Coles ou Jonathan Tawn. L’intuition est la suivante (on peut la retrouver en ligne ici). Si https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-2.2.php.png et https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-3.2.php.png ont la même loi et que l’on suppose les variables indépendantes, alors


En revanche, si les variables sont comonotones (c’est à dire égales comme on suppose les lois identiques),


Aussi, on peut supposer qu’il existe un indice https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-6.2.php.png tel que


Le soucis est que le cas d’indépendance correspond à https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-6.2.php.png=2, alors que le cas de dépendance forte correspond au cas https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-6.2.php.png=1. Il est alors usuel de faire une transformation affine pour se ramener sur [0,1], et que la force de la dépendance soit croissante avec https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-6.2.php.png, e.g.


Posons alors


qui pourra être interprété comme un (weak) tail dependence index.
Bref, ces deux mesures donnent de l’information sur le comportement dans les queues de distribution.

  • Les fonctions de concentration dans les queues

L’idée est de noter qu’il est possible d’étudier ces fonctions afin de mieux comprendre le comportement dans les queues. En s’inspirant de Gary Venter, on peut définir


pour étudier le comportement dans la queue inférieure, et


pour la queue supérieure,où https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-12.2.php.png est la copule de survie associée à https://perso.univ-rennes1.fr/arthur.charpentier/latex/toclatex2png-13.2.php.png, au sens où



Cet outil permettra de modéliser la dépendance forte. On peut également poser, afin d’étudier la dépendance faible,



  • Application statistique

L’idée est de noter qu’il est facile d’estimer ces fonctions. Ces outils peuvent être utiles pour mieux comprendre le comportement dans les queues.
Par exemple pour une copule Gaussienne de corrélation 0,5, on a la forme théorique suivante pour les fonctions de concentration (au sens fort)

Statistiquement, il est possible d’estimer ces quantités en comptant simplement le nombre d’observations dans le coin inférieur gauche, ou le coin supérieur droit.  Si on dispose d’un échantillon, on peut alors regarder ce que donnent les versions

http://perso.univ-re<br /><br /> nnes1.fr/arthur.charpentier/latex/toclatex2png-18.2.php.png



Pour un échantillon de taille n=500, on obtient les intervalles de confiance à 90% de la forme suivante,

Le code R ressemble à ça

> library(evd); data(lossalae)
> cor(lossalae,method="spearman")
         Loss     ALAE
Loss 1.000000 0.451872
ALAE 0.451872 1.000000

avec le code suivant pour la version empirique,

> z=seq(0,.5,by=.001)
> U=rank(v[,1])/(nrow(v)+1)
> V=rank(v[,2])/(nrow(v)+1)
> Lemp=rep(NA,length(z))
> Remp=rep(NA,length(z))
> for(i in 1:length(z)){
+  Lemp[i]=sum((U<=z[i])&(V<=z[i]))/sum(U<=z[i])
+  Remp[i]=sum((U>=1-z[i])&(V>=1-z[i]))/sum(U<=z[i])
+ }

et pour la version théorique,

> Lg=(pcopula(copclayton,cbind(z,z)))/(z)
> Rg=((1-2*(1-z)+pcopula(copclayton,cbind(1-z,1-z))))/(z)
> plot(c(1-z,z),c(Lg,Rg))

De plus, on a des fonctions similaires pour la dépendance au sens faible, avec le code suivant pour la version théorique,

> Lg=log(pcopula(cop,cbind(z,z)))/log(z)
> Rg=log((1-2*(1-z)+pcopula(cop,cbind(1-z,1-z))))/log(z)
> Lg=1/Lg*2-1
> Rg=1/Rg*2-1

et celui là pour la version empirique

> z=seq(0,.5,by=.001)
> v <- lossalae
> U=rank(v[,1])/(nrow(v)+1)
> V=rank(v[,2])/(nrow(v)+1)
> Lemp=rep(NA,length(z))
> Remp=rep(NA,length(z))
> for(i in 1:length(z)){
+  Lemp[i]=log(mean((U<=z[i])&(V<=z[i])))/log(mean(U<=z[i]))
+  Remp[i]=log(mean((U>=1-z[i])&(V>=1-z[i])))/log(mean(U<=z[i]))
+ }
> Lemp=1/Lemp*2-1
> Remp=1/Remp*2-1

Bref, on peut utiliser ces fonctions sur des vrais échantillons. Considérons l’exemple classique loss-alae (où l’on couple les frais dans des sinistres assurés, et les frais payés par l’assureur). On souhaite ajuster une copule, sans trop savoir laquelle. On peut commencer par étudier la dépendance forte, et comparer avec une copule Gaussienne. La copule Gaussienne de référence possède ici le même rho de Spearman que l’échantillon dont on dispose,

> cor(lossalae,method="spearman")
         Loss     ALAE
Loss 1.000000 0.451872
ALAE 0.451872 1.000000
> library(copula)
> paramgauss=.47
> paramclayton=.9
> paramgumbel=1.45
> copgauss=normalCopula(paramgauss)
> copclayton=claytonCopula(paramclayton, dim = 2)
> copgumbel=gumbelCopula(paramgumbel, dim = 2)

On obtient ici

La courbe verte est l’intervalle de confiance (ponctuel) à 95% pour une copule Gaussienne et un échantillon de même taille. On voit qu’on modélise mal la structure de dépendance. Avec une copule duale de Clayton, on obtient

et enfin pour une copule de Gumbel,

Bref, la copule de Gumbel semble réellement bien adaptée… Si on creuse en étudiant la dépendance au sens faible, on peut valider là aussi ce modèle. En effet, si la référence est la copule Gaussienne,

ou pour une copule de Clayton,

alors qu’une copule de Gumbel donnerait