Tag Archives: rstats

La pyramide des âges, au Canada

L’autre jour, Statistique Canada publiait une jolie animation interactive permettant de visualiser la déformation de la “pyramide” des âges (évoquée par exemple surhttp://www.lapresse.ca/). Bon, j’ai des données assez proches, tirées dehttp://www.mortality.org/).

> pop=read.table(
+ "http://freakonometrics.blog.free.fr/PopulationCanada.txt",
+ header=TRUE,skip=2)
> pop$Age=as.numeric(as.character(pop$Age))
> pop$Year=as.numeric(as.character(pop$Year))
> pop=pop[is.na(pop$Age)==FALSE,]

On peut faire assez facilement des pyramides des ages, par année, avec le code suivant,

> library(plotrix)
> seuils=seq(0,110,by=10)
> pop$tranche=cut(pop$Age,seuils, right = FALSE)
> an=1955
> base=pop[pop$Year==an,]
> women=aggregate(base$Female,
+ list(base$tranche),sum)[,2]
> men=aggregate(base$Female,
+ list(base$tranche),sum)[,2]
> nom=as.character(unique(pop$tranche))
> pyramid.plot(men/sum(men)*100,
+ women/sum(women)*100,labels=nom,gap=2,
+  lxcol=c("blue","blue","purple","purple","purple",
+ "purple","red","red","red","red","red"),
+  rxcol=c("blue","blue","purple","purple","purple",
+ "purple","red","red","red","red","red"))

http://freakonometrics.blog.free.fr/public/perso6/pyramide-ages-canada2.gif

Au lieu de faire des barres horizontales, on peut cumuler, par tranche d’age, et regarder l’évolution dans le temps des proportions “moins de 20 ans” (bleu), “entre 20 et 60 ans” (mauve) et “plus de 60 ans” (rouge).

> seuils=c(0,20,60,110)
> pop$tranche2=cut(pop$Age,seuils, right = FALSE)
> YEAR=1921:2008
> totaux=matrix(NA,length(YEAR),3)
> for(i in 1:length(YEAR)){
+ base=pop[pop$Year==YEAR[i],]
+ totaux[i,]=aggregate(base$Total,
+ list(base$tranche2),sum)[,2]
+ }
> stotaux=apply(totaux,1,sum)
> X=totaux[,1]/stotaux
> Y=X+totaux[,2]/stotaux

Si on regarde bien la pyramide, on s’aperçoit qu’après guerre, on a une génération importante, qui se déplace ensuite progressivement vers le haut, les baby-boomers. Pour visualiser encore davantage cette génération, on peut l’isoler sur le graphique ci-dessous, en suivant la cohorte née entre 1945 et 1950.

Mais on reviendra un peu plus en détails très bientôt sur cette génération (peut-être plutôt sur des données françaises cette fois afin de limiter les erreurs d’interprétation). A suivre…

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)
plot(focalswhite,useRaster=FALSE)

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

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
values(striped)=0
values(striped)[(h+1):(n-h)]=(values(red)[1:(n-2*h)]==
TRUE)&(values(red)[(2*h+1):n]==TRUE)
focalsstriped = focal(striped, w=3, fun=mean)
plot(focalsstriped,useRaster=FALSE)

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)
plot(focalsred,useRaster=FALSE) 
exam = stack("C:\\Users\\exam-filled.png")
red = (exam[[1]]>150)&(max(  exam[[2]] , exam[[3]])<150)
focalsred = focal(red, w=3, fun=mean)
plot(focalsred,useRaster=FALSE)

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)
plot(focalsblack,useRaster=FALSE)
correct=locator(20)
coordinates=locator(20)
X1=c(73,115,156,199,239)
X2=c(386,428.9,471,510,554)
Y=c(601,536,470,405,341,276,210,145,79,15)
LISTX=c(rep(X1,each=10),rep(X2,each=10))
LISTY=rep(Y,10)
points(LISTX,LISTY,pch=16,col="blue")

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

CORRECTX=c(X1[c(2,4,1,3,1,1,4,5,2,2)],
X2[c(2,3,4,2,1,1,1,2,5,5)])
CORRECTY=c(Y,Y)
points(CORRECTX, CORRECTY,pch=16,col="red",cex=1.3)
UNCORRECTX=c(X1[rep(1:5,10)[-(c(2,4,1,3,1,1,4,5,2,2)
+seq(0,length=10,by=5))]],
X2[rep(1:5,10)[-(c(2,3,4,2,1,1,1,2,5,5)
+seq(0,length=10,by=5))]])
UNCORRECTY=c(rep(Y,each=4),rep(Y,each=4))

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,

ind=which(values(focalsred)>.3)
yind=750-trunc(ind/610)
xind=ind-trunc(ind/610)*610
points(xind,yind,pch=19,cex=.4,col="blue")
points(CORRECTX, CORRECTY,pch=1,
col="red",cex=1.5,lwd=1.5)

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+
+ (UNCORRECTX)]
> 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…

Short versus long papers, in academic journals

This Monday, during my talk on quantile regressions (at the Montreal R-meeting), we’ve seen how those nice graphs could be interpreted, with the evolution of the slope of the linear regression, as a function of the probability level. One illustration was on large hurricanes, from Elsner, Kossin & Jagger (2008). The other one was on birthweight, from Abrevaya (2001).

It is also to illustrate that technique to academic publication, e.g. the length of papers, over time. Actually, the data we can extract from Scopus are quite similar to the ones uses on hurricanes. For several journals, it is possible to look at the length of articles. Since Scopus is quite expensive ($60,000 per year for the campus, as far as remember, so I can imagine the penalty I might have to pay for sharing such a dataset)

base=read.table("/home/scopus.csv",
header=TRUE,sep=",")
pages=base$Page.end-base$Page.start
year=base$Year

Again, a first idea can be to look at boxplots, and regression on (nonparametric) quantiles, here for Econometrica,

boxplot(pages~as.factor(year),col="light blue")
Q=function(p=.9) as.vector(by(pages,as.factor(year),
function(x) quantile(x,p)))
u=1:16
points(u,Q(p),pch=19,col="blue")
abline(lm(Q(p)~u,weights=table(year)),lwd=2,col="blue")

Consider now (as in the slides in the previous post) a quantile regression (instead of a regression on quantiles), for instance in the Annals of Probability,

library(quantreg)
u=seq(.05,.95,by=.01)
coefstd=function(u) summary(rq(pages~year,
tau=u))$coefficients[,2]
coefest=function(u) summary(rq(pages~year,
tau=u))$coefficients[,1]
CS=Vectorize(coefstd)(u)
CE=Vectorize(coefest)(u)
k=2
plot(u,CE[k,],ylim=c(min(CE[k,]-2*CS[k,]),
max(CE[k,]+2*CS[k,])))
polygon(c(u,rev(u)),c(CE[k,]+1.96*CS[k,],
rev(CE[k,]-1.96*CS[k,])),
col="light green",border=NA)
lines(u,CE[k,],lwd=2,col="red")
abline(h=0)

We have the following slope, for the year, as a function of the probability level,

The slope is always positive, so size of papers is increasing with time, short and long papers. But the influence of time is much larger for long paper than short one: for short papers (lower decile) every year, the size keeps increasing, with one more page every three years. For long paper (upper decile), it is two more pages every three years.

If we look now at the Annals of Statistics, we have

and for the evolution of the slope of the quantile regression,

Again the impact is positive: papers are longer in 2010 than 15 years ago. But the trend is the reverse: short papers (lower decile) are much longer, almost one more page every year, with long paper increase only by one more page every two years… Initially, I want to run such a study on a much longer term, with quantile regressions and splines to see when there might have been a change, both in lower and upper tails. Unfortunately, as suggested by some colleagues, there might have been some changes in the format of the journal (columns, margins, fonts, etc). That’s a shame, because I rediscover nice short papers of 5-10 pages published 20 or 30 years ago. They are nice to read (and also potentially interesting for a post on the blog). 5 pages, that’s perfect, but 40 pages, that’s way too long. I wonder if I am the only one having this feeling, missing those short but extremely interesting papers….

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,

x=c(1,2,3)
y=c(3,7,8)
plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
library(quantreg)
abline(rq(y~x,tau=.5),col="red")

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é).

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=-1,b=3.2,col="blue")

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,

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=-1,b=3.2,col="blue",lwd=2)
for(i in seq(-2,3,by=.5)){
abline(a=i,b=3.2,col="blue",lty=2)}

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)))
H=seq(-4,6,by=.01)
D=Vectorize(d)(H)
plot(H,D,type="l")

Formellement, l’optimum est ici

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

$objective
[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,

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
abline(a=optimize(d,lower=-5,upper=5)$minimum,b=3.2,
col="blue",lwd=2)
points(x[1],y[1],cex=1.8,lwd=2,col="blue")
for(i in seq(-1,5,by=.25)){
abline(a=(y[1]-i*x[1]),b=i,col="blue",lty=2)}

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)))
H=seq(-4,6,by=.01)
D=Vectorize(d2)(H)
plot(H,D,type="l")

Là encore on peut formaliser un peu

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

$objective
[1] 1.500018

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

plot(x,y,pch=19,cex=1.5,xlim=c(0,4),ylim=c(0,10))
h=optimize(d2,lower=-5,upper=5)$minimum
abline(a=(y[1]-h*x[1]),
b=h,col="purple",lwd=2)

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…

Ex-æquo aux élections ?

La saison 1 de Mad Men se termine avec l’élection de 1960, Nixon contre Kennedy, qui s’est jouée dans un mouchoir de poche.

L’autre jour, Benoit Rittaud évoquait sur http://images.math.cnrs.fr/ la possibilité que des candidats soient ex-æquo dans une élection, et les aspects légaux (car visiblement la constitution française n’a pas vraiment prévu ce cas). Le plus intéressant aurait été – à mon avis (mais mon avis est fortement biaisé) – de calculer la probabilité qu’un tel évènement survienne. Par exemple, qu’à l’issu du premier tour, le deuxième et le troisième candidat soient ex-æquo, de telle sorte qu’on ne puisse déterminer qui ira au second tour. Justement, un des commentaire affirme “la probabilité qu’un tel événement arrive est environ l’inverse du nombre typique de voix séparant le deuxième homme (ou femme) du troisième”. N’ayant réussi à démontrer ce résultat, je me suis dit qu’on pourrait se lancer dans des simulations pour quantifier cette probabilité.

Commençons par supposer qu’il y a quatre candidats (et la possibilité de s’abstenir). On suppose que le taux d’abstention est de 20%, et que les intentions de vote pour chacun des candidats soit de 20%.

EXE=function(N=100,ns=1000000){
ExE=rep(NA,ns)
for(s in 1:ns){
M=sample(c("A","B","C","D","E"),size=N,prob=rep(.2,5),
replace=TRUE)
tb=table(M)
Ms=rev(sort(tb[-(names(tb)=="A")]))
ExE[s]=(Ms[2]==Ms[3])
}
return(mean(ExE))}

Le soucis est que faire plusieurs centaines de millions de tirages de plus de 40 millions de lois multinomiales pourrait prendre du temps (genre vraiment beaucoup). Une stratégie intuitive peut être de faire des simulations pour une population plus faible, puis d’extrapoler.

E=function(N) EXE(N)
taille=c(20,50,100,200,500,1000,2000)
proba=Vectorize(E)(taille)
plot(taille,proba,type="b",log="xy")

base=data.frame(y=proba,x=taille)
reg=lm(log(y)~I(log(x)),data=base)
s2=summary(reg)$sigma
m=exp(predict(reg,newdata=
data.frame(x=41194689)))*exp(s2^2/2)

L’extrapolation linéaire n’étant peut-être pas parfaite, on peut tenter un ajustement quadratique (comme cela a été fait sur la figure ci-dessus)

reg2=lm(log(y)~I(log(x))+I(log(x)^2),data=base)
s2=summary(reg2)$sigma
m2=exp(predict(reg2,newdata=data.frame(x=41194689)))*exp(s2^2/2)

La différence de prédiction avec les deux modèles est significative comme on le voit sur le graphique. Toutefois, avoir des ex-æquo avec une très très grande population est tellement rare qu’il faut réellement beaucoup de simulations pour avoir un ordre de grandeur valide. L’extrapolation semble être une alternative valide ici (mais je suis ouvert à toute critique, ou discussion).

Maintenant que l’on a une méthodologie (ou presque), utilisons des ordres de grandeurs réalistes. Par exemple en se basant sur l’élection présidentielle de 2002. 41.19 millions d’électeurs, 28.4% d’abstentions, et pour les votes exprimés, un candidat en tête avec 19.9% des votes, un second à 16.8% et le troisième à 16.2% (et les autres plus loin). Les intentions de votes pour les principaux candidats figurent dans la loi multinomiale ci-dessous

EXE=function(N=100,ns=1000){
ExE=rep(NA,ns)
for(s in 1:ns){
M=sample(c("abs","JC","JMLP","LJ","FB","AL","JPC"
,"NM","autres"),size=N,prob=c(.284,.1988*.716,
.1686*.716,.1618*.716,.0684*.716,.0572*.716,
.0533*.716,.0525*.716,.1714104),replace=TRUE)
tb=table(M);
Ms=rev(sort(tb[-(names(tb)%in%c("abs","autres"))]));
ExE[s]=(Ms[2]==Ms[3])
}
return(mean(ExE))}

Si on en croit l’affirmation du commentaire, la probabilité d’avoir eu des ex-æquo (entre le second et le troisième candidat) serait de l’ordre de

> 1/(41194689*(.1686*.716-.1618*.716))
[1] 4.985823e-06

Je dois avouer que les calculs n’ont rien donné (on verra tout à l’heure ce qui se passe si on change la définition d’ex-æquo). Numériquement, les calculs tournent encore (après 4 jours) sur le serveur… Par contre, si la population ne contenait que 10000 électeurs, la probabilité d’avoir des ex-æquo (entre le 2ème et le 3ème candidat) serait de l’ordre de 3 sur 10000 (avec un intervalle de confiance que je n’ai toujours pas calculé)

> EXE(10000,100000000)
[1] 3.245562e-05

A titre d’illustration, on peut regarder ce qui s’est passé, commune par commune en France (les données sont en ligne sur le site du Ministère de l’Intérieur). L’exercice est biaisé parce qu’il y a de l’hétérogénéité spatiale… mais ça permet toujours d’illustrer… Les nombres de voix obtenues, pour le second et le troisième candidat, par commune est le suivant

> election=read.table(
+ "/Users/UQAM/Documents/Pres2002Tour1.csv",sep=";")
+ election=election[-1,]
> VOIX=election[,seq(14,74,by=4)]
> DEUXIEME=TROISIEME=TOTAL=rep(NA,nrow(VOIX))
> for(i in 1:nrow(VOIX)){
+ u=rev(sort(as.numeric(as.character(
+ election[i,seq(14,74,by=4)]))))
+ TOTAL[i]=sum(u)
+ DEUXIEME[i]=u[2]
+ TROISIEME[i]=u[3]
+ }

Sur 37000 communes,

> nrow(election)
[1] 37513

ou plutôt sur les 24000 communes qui ont eu plus de 5000 votes exprimés

> sum((TOTAL>5000))
[1] 24230

il y a eu 56 cas avec de deuxième et troisième candidats ex-aequo,

> base=data.frame(TOTAL,DEUXIEME,TROISIEME,COMMMUNE=election[,3])
> selection=(TOTAL>5000)
> sousbase=base[selection,]
> sum(sousbase$DEUXIEME==sousbase$TROISIEME)
[1] 56

On a même le détail,

> sousbase[indice,]
TOTAL DEUXIEME TROISIEME                     COMMMUNE
1134   5819      699       699                        Crouy
1507   5386      717       717          Louroux-Bourbonnais
2603   6235     1151      1151                  Saint-Peray
2916   6346      935       935                        Coucy
2970   7529     1413      1413                   Haraucourt
5089   5567      717       717                      Vignats

(etc). On peut alors regarder, en fonction de la taille des communes (ici le nombre de votes exprimés) la probabilité d’avoir un ex-æquo, dans une commune donnée. On va subdivisionner en 50 classes (pour commencer), et compter à chaque fois le nombre d’ex-æquo,

> n=50
> niveau=seq(1/(2*n),1-1/(2*n),by=1/n)
> base$CLASSE=cut(TOTAL,quantile(TOTAL,seq(0,1,by=1/n)),
+ labels=niveau)
> base$EXAEQUO=(base$DEUXIEME==base$TROISIEME)*1
> m=by(base$EXAEQUO, base$CLASSE, sum)

Par exemple pour les 2% des plus grosses villes, on a eu 5 cas de communes avec des ex-aequos

------------------------------------------------------------
base$CLASSE: 0.99
[1] 5

Si on regarde plus en détails, ce sont les communes avec plus de 7500 votes exprimés,

> quantile(TOTAL,.98)
98%
7427.76
> selection=(TOTAL>quantile(TOTAL,.98))
> sousbase=base[selection,]
> sum(sousbase$DEUXIEME==sousbase$TROISIEME)
[1] 5
> indice=which(sousbase$DEUXIEME==sousbase$TROISIEME)
> sousbase[indice,]
TOTAL DEUXIEME TROISIEME       COMMMUNE CLASSE EXAEQUO
2970   7529     1413      1413     Haraucourt   0.99       1
14359  7496     1392      1392          Renac   0.99       1
25028  8104     1381      1381         Heloup   0.99       1
30014  8087     1340      1340         Barnay   0.99       1
36770  7435      764       764 Pont-sur-Yonne   0.99       1

On peut représenter ces nombres de cas d’ex-æquo et faire une régression de Poisson (car on fait du comptage ici),

> taille=length(TOTAL)/n
> plot(quantile(TOTAL,niveau),m/taille)
> sb=data.frame(niveau,exaequo=as.vector(m))
> library(splines)
> reg=glm(exaequo~bs(niveau),family=poisson,data=sb)
> u=seq(0,1,by=.001)
> lines(quantile(TOTAL,u),predict(reg,
+ newdata=data.frame(niveau=u),
+ type="response")/taille,lwd=2,col="red")

pour des facilités de lecture, on a en ordonnées des probabilités et en abscisse les tailles moyennes des villes

On peut faire la même chose avec 200 classes,

Bref, pour une commune de taille moyenne (disons moyenne supérieure, avec 5000 votes exprimés), on a 2 chances sur 1000 d’avoir un second et un troisième candidat ex-æquo. Ce qui laisserait entendre que pour 40 millions, la probabilité devrait être faible, voire très faible…

Bon, au delà du problème conceptuel des ex-æquo (ou légal, car je doute qu’il soit possible d’avoir une égalité parfaite, le Ministère de l’Intérieur ferait un recompte et trouverait forcément un résultat différent), il semble que la probabilité semble beaucoup plus faible que celle obtenue par la règle ad-hoc….

Mais peut-être serait-il réaliste de dire que l’élection n’a pas su dégager une majorité (c’est le but d’une élection il me semble) si le nombre de voix séparant deux candidats est trop faible. Par exemple, si moins de 100 voix séparent les candidats, on refait un vote ! Pour revenir sur mon exemple précédant, avec 10000 votants

> EXE100(10000,1000000)
[1] 0.012641

soit une chance sur 100 (je me suis autorisé à faire moins de simulations car l’évènement me semble moins improbable… donc il y a besoin de moins de simulations pour l’observer à l’occasion). Si on va plus loin, et que l’on fait des simulations que l’on projette, comme auparavant, pour estimer la probabilité d’avoir des ex-æquo non pas pour une ville, mais au niveau de la France entière, on est plutôt aux alentours de

> (m=exp(predict(reg,newdata=
+ data.frame(x=41194689)))*exp(s2^2/2))
1
1.505422e-12

pour un modèle linéaire, et

> (m=exp(predict(reg2,newdata=
+ data.frame(x=41194689)))*exp(s2^2/2))
1
4.074588e-98

pour un modèle quadratique (avec la confiance – relativement faible – que l’on peut accorder à une telle projection). Ce qui donne une différence énorme… Si on visualise tout ça,

lines(c(1,100000000),exp(predict(reg,newdata=
data.frame(x=c(1,100000000))))*exp
(s2^2/2),lty=2,col="light blue")

que l’on représente en bleu, en vert un ajustement quadratique

lines(10^(seq(1,8,by=.1)),exp(predict(reg,newdata=
data.frame(x=10^(seq(1,8,by=.1)))))*exp
(s2^2/2),lty=2,col="dark green")

alors que le calcul suggéré dans le commentaire donnerait la courbe rouge,

lines(c(1,100000000),1/(c(1,100000000)*
(.1686*.716-.1618*.716)),col="red",lty=2)

Peut-on parler d’évènement improbable ? A titre de comparaison, quand on parle de quelque chose d’improbable, ou qui arrive avec une probabilité infinitésimale, on pense au loto. A l’ancien loto (français), il fallait choisir 6 numéros parmi 49 (sans remise). La probabilité de gagner le gros lot était de

Pour l’Euromillions, il faut choisir 5 numéros parmi 50 numéros auxquels il faut ajouter les combinaisons possibles avec les Étoiles, soit 2 numéros à choisir parmi 9. La probabilité de gagner le gros lot est alors de

Bref, tomber sur des ex-æquo lors d’une élection nationale, en France, c’est a priori plus rare que gagner le gros lot en jouant une fois au loto… On peut donc dire qu’il ne devrait pas y avoir d’ex-æquo à l’issu du premier tour !

Ruin probability and infinite time

A couple of weeks ago, I had a discussion with a practitioner, working in some financial company, about ruin, and infinite time. And it reminded me a weird result. Well, not a weird result, but a result I found disturbing, at first, when I was a student (that I rediscovered with the eyes of someone dealing with computational issues, seeing here a difficult theoretical question). Consider a simple ruin problem. A player has wealth . Then he flips a coin: tails he has a gain of 1, heads he experiences a loss of 1. At time , his wealth is where  is associated to the th coin:  is equal to 1 with probability (tails), and -1 with probability  (heads). It is also possible to write

where  can be interpreted as the net gain of the player. In order to get a good understanding of results that can be obtained. Assume  to be given. Let denote the number of heads and  the number of tails. Then , while . Let  denote the number of paths to go from point A (wealth  at time ) to point B (wealth  at time ). Note that this is a Markovian problem, that can be modeled using Markov chains

But here, we will focus on combinatorial results. Hence,

In order to derive probabilities to reach , let  denote the number of paths going from  to . And let denote the number of paths going from  to  that do reach  at some point between  and . Using a simple reflexion property, then if  and  are positive,

Based on those reflexions, two results can be derived (focusing on probability, instead of counting paths). First, we can obtain that

(given that n and x have the same parity). The second result we can obtain is that

Based on those two expressions, if  denotes the first time  become null, given ,

then

This can be computed easily,

> x=10
> p=.55
> ProbN=function(n){
+ pb=0
+ if(abs(n-x) %% 2 == 0)
+ pb=x/n*choose(n,(n+x)/2)*(1-p)^((n+x)/2)*(p)^((n-x)/2)
+ return(pb)}
> plot(Vectorize(ProbN)(1:1000),type="s")

That looks nice… But if we look closer, we can wonder what

would be ? Since we have the distribution of a probabilty measure, we might expect one. But here

> sum(Vectorize(ProbN)(1:1000))
[1] 0.134385

And this is not due to calculation mistakes that we do not get 1 here. Actually, we should write

which might be interpreted as the probability of ruin, starting from , that we denote  from now on. The term on the left can be approximated using monte-carlo simulations

> p=.55
> x=10
> m=1000
> simul=10000
> S=sample(c(-1,1),size=m*simul,replace=TRUE,prob=c(1-p,p))
> MS=matrix(S,simul,m)
> for(k in 2:m) MS[,k]=MS[,k]+MS[,k-1]
> T0=function(vm) which(vm<=(-x))[1]
> MTmin=apply(MS,1,T0)
> mean(is.na(MTmin)==FALSE)
[1] 0.1328

To check the validity of the relationship above, a simple (theoretical) recursive formula can be derived for the term on the right (ruin probability), namely

with a boundary conditions , and . Then is comes that

Note that it might be tricky to check using monte carlo simulation… since we cannot have an infinite number of runs. And we’re dealing precisely with things that do occur when time is infinite. Actually, we can still check convergence, considering an upper limit  for the number of runs, and then letting  go to infinity. Note that an explicit formula can then be derived (using additional border condition )

Using the following code, it is possible to calculate ruin probability, in order to estimate .

> MSmin=apply(MS,1,min)
> mean(MSmin<=(-x))
[1] 0.1328
> (((1-p)/p)^x-((1-p)/p)^m)/(1-((1-p)/p)^m)
[1] 0.1344306

The following graph shows the evolution of ruin probability as a function of initial wealth (with monte carlo simulation, with a fixed horizon – including a confidence interval – versus the analytical expression)

Hence, with stopping times, one should remember that

and that those two terms can be approximated simply using simulations or standard approximations.

Simulation d’un processus de Lévy, et discrétisation

Avec @renaudjf, on discutait l’autre jour de la simulation d’un processus de Lévy. Et on se posait la question d’un algorithme optimal pour combiner un processus de Poisson (ou un process Poisson composé) avec un processus de Wiener (avec éventuellement un drift, voire une diffusion plus générale). En fait, pour générer des processus de Poisson, j’ai toujours eu l’habitude de simuler les durées entre sauts (avec des lois exponentielles, indépendantes, comme dans des vieux billets). Jean François me suggérait d’utiliser une propriété d’uniformité des sauts sur un intervalle de temps donné, conditionnellement aux nombres de sauts.

Commençons par la première piste. On peut générer un processus de Wiener, éventuellement avec un drift, et à coté, on peut générer les lois exponentielles  (qui vont correspondre aux durées entre sauts), et éventuellement aussi des amplitudes de sauts (e.g. des pertes qui suivent des lois exponentielles). On a ici

où . On commence par générer  en notant que

où les incréments  sont Gaussiens (centrés et de variance ) et indépendants les uns des autres. Quant aux durées entre sauts, les  , ce sont des lois exponentielles indépendantes, de moyenne . Voilà le code qui permet de générer les trois composantes,

n=1000
h=1/n
lambda=5
set.seed(2)
W=c(0,cumsum(rnorm(n,sd=sqrt(h))))
W=rexp(100,lambda)
N=sum(cumsum(W)<1)
T=cumsum(W[1:N])
X=-rexp(N)

Le hic est que pour le processus de Wiener, on a du discrétiser, alors que pour le processus de Poisson composé, non. Pourtant, il va bien se ramener sur une même échelle de temps. Une première piste est de créer vraiment la fonction 

Lt=function(t){
W[trunc(n*t)+1]+sum(X[T<=t])+lambda*t
}

et pour faire un dessin ensuite, c’est un jeu d’enfant

L=Vectorize(Lt)
u=seq(0,1,length=n+1)
plot(u,L(u),type="l",col="blue")

Une autre possibilité est d’utiliser une propriété d’uniformité du processus de Poisson que j’évoquais en introduction. Car le processus de Poisson vérifie une propriété remarquable: si  est la date où survient le ième saut, , alors conditionnellement au fait que , les variables  correspondent aux statistiques d’ordres de  variables indépendantes, uniformément distribuées sur , i.e.

Cette propriété se trouve dans Wolff (1982). L’idée de la démonstration est relativement simple. Commençons avec un (unique) saut, alors pour ,

i.e. on retrouve la fonction de répartition d’une loi uniforme sur . On itère ensuite avec 2 sauts, 3 sauts, etc.

La traduction en R de cette idée est tout simplement (car on se place sur )

N=rpois(1,lambda)
T=runif(N)
X=-rexp(N)

Ensuite, une stratégie est de discrétiser le processus de Poisson, avec le même pas de temps que le processus de Wiener,

indice=trunc(T*n)+1
saut=rep(0,n+1)
saut[indice]=X
processus=W+cumsum(saut)+lambda*u

On retrouve la même trajectoire qu’auparavant

plot(u,processus,type="l",col="red")

,,

Sauf qu’on a eu de la chance. Avec cette procédure, il ne faut pas que l’on ait deux sauts dans le même intervalle de temps ! Bon, il est vrai qu’une caractérisation du processus de Poisson est que

et donc on doit avoir très peu de chance d’avoir deux sauts au même instant d’autant plus que le pas de temps est petit. Mais “peu de chance” ne veut pas dire nul, et si on génère des milliers de trajectoires, la probabilité d’avoir une fois un soucis n’est pas négligeable.

Jean-François a eu l’idée brillante de proposer de tirer non pas des lois uniformes sur , mais des lois uniformes discrètes, dans

T=c(0,sort(sample((1:(n-1)/n),size=N,replace=FALSE)))

sans remise afin d’éviter d’avoir deux sauts au même moment. L’idée était séduisante, mais il me restait à etre convaincu que les durées entre sauts… continuaient à être (presque) exponentielles.

Pour ça, on peut faire quelques tests. Par exemple, générer quelques simulations pour avoir une centaines de sauts (et donc une centaine de durées entre sauts), puis faire un test de loi exponentielle (de moyenne )

VT=0
for(ns in 1:20){
N=rpois(1,lambda)
if(N>0){
T=c(0,sort(sample((1:(n-1)/n),size=N,replace=FALSE)))
VT=c(VT,diff(T))
}}

On fait ici 20 boucles car on avait fixé

lambda=5

et j’avais dit que je voulais une centaine d’observations pour faire un test de loi (ce qui est purement arbitraire, on en conviendra). On peut ensuite faire un test d’ajustement de loi exponentielle,

ks.test(VT[-1],"pexp",lambda)$p.value

Si on répète un grand nombre de fois, en changeant le pas de temps (ou le nombre de subdivisons de l’intervalle de temps), on notera qu’effectivement, avec un grand pas de temps (à gauche ci-dessous) on rejette souvent, voire presque toujours, l’hypothèse de loi exponentielle. Mais qu’assez vite, c’est une hypothèse qui n’est pas invraisemblable,

Je ne sais pas entre ces deux fonctions ce qui serait le plus rapide, mais on a deux jolis algorithmes pour générer les processus de Lévy. Non ?

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.

http://freakonometrics.hypotheses.org/files/2015/12/lex03.gif

http://freakonometrics.hypotheses.org/files/2015/12/le04.gif désigne un poids attribué à la nouvelle observation dans la fonction de lissage http://freakonometrics.hypotheses.org/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")

http://freakonometrics.hypotheses.org/files/2015/12/lissage-exp-1.gif

Sur la figure ci-dessus, non visualise l’impact de http://freakonometrics.hypotheses.org/files/2015/12/lisse16.gif sur le lissage. La prévision que l’on fait à la date http://freakonometrics.hypotheses.org/files/2015/12/lex14.gif, pour un horizon http://freakonometrics.hypotheses.org/files/2015/12/lisse12.gif est alors http://freakonometrics.hypotheses.org/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

http://freakonometrics.hypotheses.org/files/2015/12/lisse14.gif

> 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.

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

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

Coefficients:
[,1]
a 805.0389

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

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

http://freakonometrics.hypotheses.org/files/2015/12/lex50.gif

et

http://freakonometrics.hypotheses.org/files/2015/12/HW02.gif

Dans ce cas, la prédiction est http://freakonometrics.hypotheses.org/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 http://freakonometrics.hypotheses.org/files/2015/12/lisse15.gif (http://freakonometrics.hypotheses.org/files/2015/12/lisse16.gif étant ici fixé),
http://freakonometrics.hypotheses.org/files/2015/12/lissage-exp-2.gif
(le lissage simple – avec le même poids http://freakonometrics.hypotheses.org/files/2015/12/lisse16.gif – apparaissant ici en trait clair).

Fisher-Tippett theorem and limiting distribution for the maximum

Tomorrow, we will discuss Fisher-Tippett theorem. The idea is that there are only three possible limiting distributions for normalized versions of the maxima of i.i.d. samples https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-00.gif. For bounded distribution, consider e.g. the uniform distribution on the unit interval, i.e. https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-09.gif on the unit interval. Let https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-10.gif and https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-11.gif. Then, for all https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-12.gif and https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-13.gif,

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-14.gif

i.e. the limiting distribution of the maximum is Weibull’s.

set.seed(1)
s=1000000
n=100
M=matrix(runif(s),n,s/n)
V=apply(M,2,max)
bn=1
an=1/n
U=(V-bn)/an
hist(U,probability=TRUE,,col="light green",
xlim=c(-7,1),main="",breaks=seq(-20,10,by=.25))
u=seq(-10,0,by=.1)
v=exp(u)
lines(u,v,lwd=3,col="red")

For heavy tailed distribution, or Pareto-type tails, consider Pareto samples, with distribution function https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-05.gif. Let https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-06.gif and https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-07.gif, then

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-08.gif

which means that the limiting distribution is Fréchet’s.

set.seed(1)
s=1000000
n=100
M=matrix((runif(s))^(-1/2),n,s/n)
V=apply(M,2,max)
bn=0
an=n^(1/2)
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(0,7),main="",breaks=seq(0,max(U)+1,by=.25))
u=seq(0,10,by=.1)
v=dfrechet(u,shape=2)
lines(u,v,lwd=3,col="red")

For light tailed distribution, or exponential tails, consider e.g. a sample of exponentially distribution variates, with common distribution function https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-01.gif. Let https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-02.gif and https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-03.gif, then

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-04.gif

i.e. the limiting distribution for the maximum is Gumbel’s distribution.

library(evd)
set.seed(1)
s=1000000
n=100
M=matrix(rexp(s,1),n,s/n)
V=apply(M,2,max)
(bn=qexp(1-1/n))
log(n)
an=1
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(-2,7),ylim=c(0,.39),main="",breaks=seq(-5,15,by=.25))
u=seq(-5,15,by=.1)
v=dgumbel(u)
lines(u,v,lwd=3,col="red")

Consider now a Gaussian https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-17.gif sample. We can use the following approximation of the cumulative distribution function (based on l’Hopital’s rule)

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-15.gif

as https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-16.gif. Let https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-18.gif and https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-19.gif. Then we can get

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-20.gif

as https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-21.gif. I.e. the limiting distribution of the maximum of a Gaussian sample is Gumbel’s. But what we do not see here is that for a Gaussian sample, the convergence is extremely slow, i.e., with 100 observations, we are still far away from Gumbel distribution,

and it is only slightly better with 1,000 observations,

set.seed(1)
s=10000000
n=1000
M=matrix(rnorm(s,0,1),n,s/n)
V=apply(M,2,max)
(bn=qnorm(1-1/n,0,1))
an=1/bn
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(-2,7),ylim=c(0,.39),main="",breaks=seq(-5,15,by=.25))
u=seq(-5,15,by=.1)
v=dgumbel(u)
lines(u,v,lwd=3,col="red")

Even worst, consider lognormal observations. In that case, recall that if we consider (increasing) transformation of variates, we are in the same domain of attraction. Hence, since https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-22.gif, if

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-23.gif

then

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-24.gif

i.e. using Taylor’s approximation on the right term,

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2018/02/max-25.gif

This gives us normalizing coefficients we should use here.

set.seed(1)
s=10000000
n=1000
M=matrix(rlnorm(s,0,1),n,s/n)
V=apply(M,2,max)
bn=exp(qnorm(1-1/n,0,1))
an=exp(qnorm(1-1/n,0,1))/(qnorm(1-1/n,0,1))
U=(V-bn)/an
hist(U,probability=TRUE,col="light green",
xlim=c(-2,7),ylim=c(0,.39),main="",breaks=seq(-5,40,by=.25))
u=seq(-5,15,by=.1)
v=dgumbel(u)
lines(u,v,lwd=3,col="red")

Construire une courbe ROC

Juste avant les vacances, Jean-Pierre Liégeois, un jeune lecteur du var, me demandais par courriel, “à partir d’une régression logistique (ou d’une matrice de confusion 2×2), comment programmer en R, un programme qui construit la courbe ROC associée“. Avant d’aller plus loin (et de répondre a la question), je vais renvoyer vers un vieux billet sur les matrices de confusion. L’idée est que l’on suppose que l’on dispose d’un prédicteur d’une variable prenant des valeurs 0 et 1 (ou pour reprendre la terminologie classique “positif” et “négatif”), par exemple un modèle logistique. Formellement, pour l’ensemble de nos observations, on a une valeur observée http://freakonometrics.hypotheses.org/files/2018/02/ROC-01.png et (comme je l’expliquais dans un autre billet) et d’un score \widehat{S}. Et c’est ce score qu’on va utiliser pour construire la courbe ROC. Ce score sera utilise pour prédire http://freakonometrics.hypotheses.org/files/2018/02/ROC-02.png . La règle d’affectation est alors simple: on se fixe un seuil http://freakonometrics.hypotheses.org/files/2018/02/ROC-04.png, et

  • si http://freakonometrics.hypotheses.org/files/2018/02/ROC-05.png, alors http://freakonometrics.hypotheses.org/files/2018/02/ROC-02.png est “positif”
  • si http://freakonometrics.hypotheses.org/files/2018/02/ROC-06.png, alors http://freakonometrics.hypotheses.org/files/2018/02/ROC-02.png est “négatif”

On peut alors construire une matrice dite de confusion, qui est simplement un table de contingence,

valeur observée http://freakonometrics.hypotheses.org/files/2018/02/ROC-01.png
valeur prédite
http://freakonometrics.hypotheses.org/files/2018/02/ROC-02.png
“positif” “négatif”
“positif” TP FP
“négatif” FN TN

où TP désigne les vrais positifs (true positive), TN les vrais négatifs (true negative),FP désigne les faux positifs (false positive) ou erreur de type I (dans une terminologie de théorie de la décision, ou de théorie des tests), et FN désigne les faux négatifs (false negative) ou erreur de type II.
Quid de la mise en œuvre sous R ? Commençons par générer des données, et estimons un modèle de régression.

set.seed(1)
n=50
X=rnorm(n)
Y=rbinom(n,size=1,prob=
exp(2*X-1)/(1+exp(2*X-1)))
B=data.frame(Y,X)
reg=glm(Y~X,family=binomial,data=B)
S=predict(reg,type="response")

On a maintenant notre observation (variable prenant les valeurs 0 ou 1) et notre score. On va ensuite pouvoir choisir plusieurs valeurs possible pour le seuil, et visualiser le taux de vrais positifs, en fonction du taux de faux positifs.

plot(0:1,0:1,xlab="False Positive Rate",
ylab="True Positive Rate",cex=.5)
for(s in seq(0,1,by=.01)){
Ps=(S>s)*1
FP=sum((Ps==1)*(Y==0))/sum(Y==0)
TP=sum((Ps==1)*(Y==1))/sum(Y==1)
points(FP,TP,cex=.5,col="red")
}

On a alors le graphique suivant,

Si on relit les points, on a alors la courbe ROC,

FP=TP=rep(NA,101)
plot(0:1,0:1,xlab="False Positive Rate",
ylab="True Positive Rate",cex=.5)
for(s in seq(0,1,by=.01)){
Ps=(S>s)*1
FP[1+s*100]=sum((Ps==1)*(Y==0))/sum(Y==0)
TP[1+s*100]=sum((Ps==1)*(Y==1))/sum(Y==1)
}
lines(c(FP),c(TP),type="s",col="red")

En fait, le code est assez simple, et il traîne dans différents packages de R, e.g.

library(ROCR)
pred=prediction(P,Y)
perf=performance(pred,"tpr", "fpr")
plot(perf,colorize = TRUE)

On peut aussi s’amuser a bootstrapper l’échantillon pour construire des intervalles de confiance, ou ajuster des modèles théoriques,

library(verification)
roc.plot(Y,P, xlab = "False Positive Rate",
ylab = "True Positive Rate", main = "", CI = TRUE,
n.boot = 100, plot = "both", binormal = TRUE)

ou encore (toujours avec des bornes de confiance obtenues par bootstrap)

library(pROC)
PROC=plot.roc(Y,P,main="", percent=TRUE,
ci=TRUE)
SE=ci.se(PROC,specificities=seq(0, 100, 5))
plot(SE, type="shape", col="light blue")

Maps with R, and polygon boundaries

With R, it is extremely easy to draw maps. Let us start with something simple, like French regions. Baptiste mentioned on his blog that shapefiles can be downloaded from http://ign.fr/ website. Hence, if you extract the zip file, it is possible to get claims frequency per region (as done in the course ACT2040),

> library(maptools)
> library(maps)
> departements<-readShapeSpatial("DEPARTEMENT.SHP")
> region<-tapply(baseFREQ[,"nbre"],
+ as.factor(baseFREQ[,"region"]),sum)/
+ tapply(baseFREQ[,"exposition"],
+ as.factor(baseFREQ[,"region"]),sum)
> depFREQ=rep(NA,nrow(departements))
> names(depFREQ)=as.character(
+ departements$CODE_REG)
> for(nom in names(region)){
+ depFREQ[names(depFREQ)==nom] =
+ region[nom]}
> plot(departements,col=gray((depFREQ-.05)*20))
> legend(166963,6561753,legend=seq(1,0,by=-.1)/20+.05,
+ fill=gray(seq(1,0,by=-.1)),cex=1.25, bty="n")

Another application is on earthquakes. It is possible to use shapefiles of tectonic plates contour, and to relate earthquakes to plates. Shapefiles can be found onhttp://www.colorado.edu/ (here).

http://freakonometrics.blog.free.fr/public/perso4/plate-tekto.gif

First, we can extract the shapes of the tectonic plates

> plates = readShapePoly("plates.shp",
+ proj4string=CRS("+proj=longlat"))
> PP=SpatialPolygons2PolySet(plates)

Consider Montreal,

> montreal=c(-73.600,45.500)

Given that specific location, it is possible to use the following code to get the associated plate,

> PLATE.loc=function(pt){
+ K=NA
+ for(k in 1:17){
+ c=point.in.polygon(pt[1], pt[2],
+ PP[PP$PID==k,c("X")],PP[PP$PID==k,c("Y")],
+ mode.checked=FALSE)
+ if(c>0){K=k}
+ }
+ return(K)}
> abline(v=montreal[1],col="red")
> abline(h=montreal[2],col="red")
> PLATE.loc(montreal)
[1] 1

and then to plot the associated tectonic plate very easily

> PLATE=function(k0){
+ library(maps)
+ map("world")
+ polygon(PP[PP$PID==k0,c("X")],PP[PP$PID==k0,c("Y")],
+ col="red")
+ for(k in (1:17)[-k0]){polygon(PP[PP$PID==k,c("X")],
+ PP[PP$PID==k,c("Y")],col="light blue")}
+ map("world",add=TRUE)}
> PLATE(PLATE.loc(montreal))

Those code were used in the paper written with Mathieu, and that will be presented on January 30th at the Geotop seminar.

Basics on Markov Chain (for parents)

Markov chains is a very interesting and powerful tool. Especially for parents. Because if you think about it quickly, most of the games our kids are playing at are Markovian. For instance, snakes and ladders…

It is extremely easy to write down the transition matrix, one just need to define all snakes and ladders. For the one above, we have,

n=100
M=matrix(0,n+1,n+1+6)
rownames(M)=0:n
colnames(M)=0:(n+6)
for(i in 1:6){diag(M[,(i+1):(i+1+n)])=1/6}
M[,n+1]=apply(M[,(n+1):(n+1+6)],1,sum)
M=M[,1:(n+1)]
starting=c(4,9,17,20,28,40,51,54,62,
64,63,71,93,95,92)
ending  =c(14,31,7,38,84,59,67,34,19,
60,81,91,73,75,78)
for(i in 1:length(starting)){
v=M[,starting[i]+1]
ind=which(v>0)
M[ind,starting[i]+1]=0
M[ind,ending[i]+1]=M[ind,ending[i]+1]+v[ind]}

So, why is it important to have a Markov Chain ? Because, once you’ve noticed that you had a Markov Chain game, you can derive anything you want. For instance, you can get the distribution after some turns,

powermat=function(P,h){
Ph=P
if(h>1){
for(k in 2:h){
Ph=Ph%*%P}}
return(Ph)}
initial=c(1,rep(0,n))
COLOR=rev(heat.colors(101))
u=1:sqrt(n)
boxes=data.frame(
index=1:n,
ord=rep(u,each=sqrt(n)),
abs=rep(c(u,rev(u)),sqrt(n)/2))
position=function(h=1){
D=initial%*%powermat(M,h)
plot(0:10,0:10,col="white",axes=FALSE,
xlab="",ylab="",main=paste("Position after",h,"turns"))
segments(0:10,rep(0,11),0:10,rep(10,11))
segments(rep(0,11),0:10,rep(10,11),0:10)
for(i in 1:n){
polygon(boxes$abs[i]-c(0,0,1,1),
boxes$ord[i]-c(0,1,1,0),
col=COLOR[min(1+trunc(500*D[i+1]),101)],
border=NA)}
text(boxes$abs-.5,boxes$ord-.5,
boxes$index,cex=.7)
segments(c(0,10),rep(0,2),c(0,10),rep(10,2))
segments(rep(0,2),c(0,10),rep(10,2),c(0,10))}

Here, we have the following (note that I assume that once 100 is reached, the game is over)

Assume for instance, that after 10 turns, your daughter accidentally drops her pawn out of the game. Here is the theoretical (unconditional) position of her pawn after 10 turns,

 so, if she claims she was either on 58, 59 or 60, here are the theoretical probabilities to be in each cell after 10 turns,

> h=10
> (initial%*%powermat(M,h))[59:61]/
+ sum((initial%*%powermat(M,h))[59:61])
[1] 0.1597003 0.5168209 0.3234788

i.e. it is more likely she was on 59 (60th cell of the vector since we start in 0). You can also look at the distribution of the number of turns (at first with only one player).

distrib=initial%*%M
game=rep(NA,1000)
for(h in 1:length(game)){
game[h]=distrib[n+1]
distrib=distrib%*%M}
plot(1-game[1:200],type="l",lwd=2,col="red",
ylab="Probability to be still playing")

Once you have that survival distribution, you have the expected number of turns to finish the game,

> sum(1-game)
[1] 32.16499

i.e. in 33 turns, on average, your daughter reaches the 100 cell. But in 50% of the games, it takes less than 29,

> max(which(1-game>.5))
[1] 29

But assuming that you are playing with your daughter, and that the game is over once one player reaches the 100 cell, it is possible to get the survival distribution of the first time one of us reaches the 100 cell,

plot((1-game[1:200])^2,type="l",lwd=2,col="blue",
ylab="Probability to be still playing (2 players)")

Here, the expected number of turns before ending the game is

> sum((1-game)^2)
[1] 23.40439

And if you ask your son to join the game, the survival distribution function is

plot((1-game[1:200])^3,type="l",lwd=2,col="purple",
ylab="Probability to be still playing (3 players)")

i.e. the expected number of turns before the end is now

> sum((1-game)^3)
[1] 20.02098

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

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

puis 1960-2000

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

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

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

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

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

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

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

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

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

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

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

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

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

ACT2040: examen final (suite et fin)

Comme annoncé, l’examen terminal a lieu ce matin, et il est basé sur des données mises en ligne sur le blog depuis une dizaine de jours. L’énoncé est un peu longue, avec une bonne trentaine de pages de sorties informatiques, et un énoncécomprenant une trentaine de questions. Des éléments de correction sont en ligne ci-dessous (plutôt que de taper une correction en pdf, autant mettre le code en ligne).

On suppose importée les bases (sinon, le code est ici). Pour la première partie – sur la tarification (a priori) – la fréquence annuelle de sinistres est tout simplement

> mean(BASEN$nombre)/mean(BASEN$exposition)
[1] 0.05241692
> exp(glm(nombre~1,offset=log(exposition),data=BASEN,
+ family=poisson(link="log"))$coefficients)
(Intercept)
0.05241692

(je vais essayer de donner à chaque question deux méthodes pour obtenir le résultat). Si on suppose que le nombre de sinistres par an suit une loi de Poisson, alors la probabilité d’avoir au moins un accident est

> lambda=mean(BASEN$nombre)/mean(BASEN$exposition)
> 1-ppois(0,lambda)
[1] 0.05106684
> 1-exp(-lambda)
[1] 0.05106684

Si on cherche à distinguer par sexe du conducteur principal

> mean(BASEN$nombre[BASEN$sexeconducteur=="H"])/
+ mean(BASEN$exposition[BASEN$sexeconducteur=="H"])
[1] 0.05218143
> mean(BASEN$nombre[BASEN$sexeconducteur=="F"])/
+ mean(BASEN$exposition[BASEN$sexeconducteur=="F"])
[1] 0.05457013
> exp(glm(nombre~0+sexeconducteur,offset=log(exposition),
+ data=BASEN,family=poisson(link="log"))$coefficients)
sexeconducteurF sexeconducteurH
0.05457013      0.05218143

Ces deux valeurs ne sont pas significativement différentes. Si on regarde la sortie de la régression

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
sexeconducteurF -2.90827    0.05464  -53.23   <2e-16 ***
sexeconducteurH -2.95303    0.01848 -159.82   <2e-16 ***

les deux coefficients sont quasiment identiques (en enlevant la constante). Si on rajoute la constante, la modalité de référence devientle sexe féminin

Coefficients:
                Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -2.90827    0.05464 -53.230   <2e-16 ***
sexeconducteurH -0.04476    0.05768  -0.776    0.438

et la modalité homme n’est alors plus significative. Bref, le sexe du conducteur n’est pas une variable explicative.

Ensuite, parmi les variables qui semble pouvoir etre utilisée pour modéliser le nombre de sinistre, on pourra retenir l’age du conducteur (décroissant jusqu’à 30 ans, puis croissant), l’age du véhicule (décroissant), l’ancienneté du permis de conduire (décroissant les 10 premières années puis croissant – ou indépendant), éventuellement certains type de véhicules (le coupé cabriolet semble avoir – en moyenne – plus d’accident que la Jeep), la situation familiale, le lieu d’habitation (en particulier rural semble avoir moins d’accident qu’urbain), le statut d’habitation, le type de paiement de la prime (mensuel semblant plus risqué qu’annuel), le poids du véhicule, ou la marque (ou plutot certaines marques, comme Peugeot versus Burstner Mobil). Bref, toutes, sauf le sexe du conducteur.

Pour un conducteur de 45 ans qui a son permis de conduire depuis seulement 5 ans, les trois modèles donnent la fréquence annuelle suivante

> nbase=data.frame(ageconducteur=45,agepermis=5,exposition=1)
> exp(t(c(1,5))%*%reg6$coefficients)
[,1]
[1,] 0.03395972
> predict(reg6,newdata=nbase,type="response")
1
0.03395972
> exp(t(c(1,45))%*%reg7$coefficients)
[,1]
[1,] 0.04100356
> predict(reg7,newdata=nbase,type="response")
1
0.04100356
> exp(t(c(1,5,45))%*%reg8$coefficients)
[,1]
[1,] 0.05395398
> predict(reg8,newdata=nbase,type="response")
1
0.05395398

Les trois modèles sont respectivement un modèle de Poison avec un lien log, un modèle binomiale négative avec un lien log, et un modèle binomial, avec un lien logit. Les deux premiers sont très proches: ils ont la meme fonction de lien, et comme nous l’avions vu à l’examen intra, le modèle binomiale négative est un modèle de Poisson avec une variable latente non observée qui suivrait une loi Gamma. Donc pour ces deux modèles, si les classes sont homogènes, et que par classe, la loi de Poisson est adaptée, il est normal d’avoir des coefficients très proches. Pour le modèle binomial avec un lien logistique, il faut se souvenir que la fonction lien est

http://freakonometrics.blog.free.fr/public/perso4/act2040-1.gif

car ici la probabilité d’avoir un accident est faible. Donc les deux fonctions liens sont proches. De plus, pour une loi de Poisson

http://freakonometrics.blog.free.fr/public/perso4/act2040-2.gif
http://freakonometrics.blog.free.fr/public/perso4/act2040-3.gif

Bref, pour notre loi de Poisson, on a presque une loi binomiale de paramètre http://freakonometrics.blog.free.fr/public/perso4/act2040-4.gif si http://freakonometrics.blog.free.fr/public/perso4/act2040-4.gif est petit.

Pour la prédiction suivant les trois modèles, on utilise classiquement

> exp(t(c(1,40,0,0,1,0,1,0))%*%reg9$coefficients)
           [,1]
[1,] 0.04351037

Comme on a accès aux sorties, pour les trois modèles, les trois prédictions sont très proches (mais on s’en doutait d’après la question précédante)

> predict(reg9,newdata=nbase,type="response")
         1 
0.04351037 
> predict(reg10,newdata=nbase,type="response")
         1 
0.04351056 
> predict(reg11,newdata=nbase,type="response")
         1 
0.04352741

Pour la méthode de biais minimal, il s’agit du principe de balancement, qui coïncide avec la régression de Poisson sur les classes. Pour les régressions stepwise, dans un cas, on élimine (étape par étape) des variables, et dans l’autre cas, des modalités de variables (pour les variables qualitatives).

Les deux modèles semblent valides, au sens où toutes les variables sont significatives. Le découpage en classes pour l’age du véhicule semble justifié par l’arbre de régression, qui suggèrait de découper à 8 et 14 ans. Pour les prédictions, on a pour le premier modèle

> nbase15=data.frame(ageconducteur=40,agepermis=20,
+ habitation="urbain",locataire=TRUE,paymensuel=FALSE,
+ marie=TRUEFORD=TRUE,RENAULT=FALSE,PEUGEOT=FALSE,
+ exposition=1)
> predict(reg15,newdata=nbase15,type="response")
         1 
0.03398312

alors que le second donne

> nbase16=data.frame(ageconducteur=40,agevehF="NEUF",rural=FALSE,
+ locataire=TRUE,paymensuel=FALSE,marie=TRUE,
+ jeune=FALSE,ancien=FALSE,exposition=1)
> predict(reg16,newdata=nbase16,type="response")
         1 
0.05208656

Moralité, nos deux modèles (qui semblent autant valides l’un que l’autre) donnent des prédictions relativement différentes. Si l’assuré change de voiture, cela n’a pas d’influence pour le second modèle (le véhicule est toujours dans la catégorie neuf). Pour le premier, il perd le rabais qu’il avait en conduisant un véhicule Ford, et donc la fréquence annuelle de sinistre augmente de 23.5%

> 1/exp(reg15$coefficients["FORDTRUE"])
FORDTRUE 
1.235329

Bref, la prédiction devient

> predict(reg15,newdata=nbase15,type="response")/
+ exp(reg15$coefficients["FORDTRUE"])
         1 
0.04198033

On peut d’ailleurs le vérifier simplement,

> nbase15b=data.frame(ageconducteur=40,agepermis=20,
+ habitation="urbain",locataire=TRUE,paymensuel=FALSE,
+ marie=TRUE,FORD=FALSE,RENAULT=FALSE,PEUGEOT=FALSE,
+ exposition=1)
> predict(reg15,newdata=nbase15b,type="response")
         1 
0.04198033

Si on passe maintenant à la modélisation de la charge de sinistres, le coût moyen par sinistre serait

> mean(BASEY$cout)
[1] 1468.522

Cette fois, contrairement à la modélisation de la fréquence, il est plus difficile de trouver des variables explicatives. A priori, au vu des graphiques, aucune variable ne semblerait convenir…
Si on regarde la prédiction pour le premier modèle obtenu,

> nbase7=data.frame(VOLKSWAGEN=FALSE,agepermis=20,
+ sexeconducteur="H",usage="TOUS_DEPLACEMENTS")
> predict(regr7,newdata=nbase7,type="response")
       1 
1388.622

alors que pour le second,

> exp(predict(regr8,newdata=nbase7))*
+ exp(summary(regr8)$sigma^2/2)
       1 
1359.552

Si c’est la femme de l’assuré qui avait pris l’assurance, cela n’aurait rien changé pour le second modèle. En revanche, pour le premier, l’assuré perdrait le rabais qu’il avait en étant un homme, et donc le cout moyen attendu augmenterait de 23%,

> 1/exp(regr7$coefficients["sexeconducteurH"])
sexeconducteurH 
       1.237832

i.e.

> predict(regr7,newdata=nbase7,type="response")/
+ exp(regr7$coefficients["sexeconducteurH"])
       1 
1718.88

La encore, on retrouve le meme résultat directement

> nbase7b=data.frame(VOLKSWAGEN=FALSE,agepermis=20,
+ sexeconducteur="F",usage="TOUS_DEPLACEMENTS")
> predict(regr7,newdata=nbase7b,type="response")
       1 
1718.88

Passons maintenant aux calculs de la prime pure. Sans segmentation, la prime pure serait

> sum(BASEN$nombre)/sum(BASEN$exposition)*
+ mean(BASEY$cout)
[1] 76.9754

Les prédictions des modèles pour la fréquence donnaient

> (n1=predict(reg15,newdata=nbase15,type="response"))
         1 
0.03398312 
> (n2=predict(reg16,newdata=nbase16,type="response"))
         1 
0.05208656

alors que pour les couts moyens par sinistre, nous avions

> (c1=predict(regr7,newdata=nbase7,type="response"))
       1 
1388.622 
> (c2=exp(predict(regr8,newdata=nbase7))
+ *exp(summary(regr8)$sigma^2/2))
       1 
1359.552

Les quatre primes pures obtenues sont alors

> c(n1,n2)%*%t(c(c1,c2))
            1        1
[1,] 47.18972 46.20183
[2,] 72.32855 70.81440

Passons maintenant à la seconde partie – sur le provisionnement des sinistres à payer. On commence par travailler sur le triangle

TRIANGLE2010
          D1      D2      D3      D4      D5      D6      D7      D8      D9
A2002 357848 1124788 1735330 2218270 2745596 3319994 3466336 3606286 3833515
A2003 352118 1236139 2170033 3353322 3799067 4120063 4647867 4914039 5339085
A2004 290507 1292306 2218525 3235179 3985995 4132918 4628910 4909315      NA
A2005 310608 1418858 2195047 3757447 4029929 4381982 4588268      NA      NA
A2006 443160 1136350 2128333 2897821 3402672 3873311      NA      NA      NA
A2007 396132 1333217 2180715 2985752 3691712      NA      NA      NA      NA
A2008 440832 1288463 2419861 3483130      NA      NA      NA      NA      NA
A2009 359480 1421128 2864498      NA      NA      NA      NA      NA      NA
A2010 376686 1363294      NA      NA      NA      NA      NA      NA      NA

Pour obtenir le montant payé en 2008 pour les sinistres survenus en 2004, on utilise la différence entre

> TRIANGLE2010["A2004",c("D4","D5")]
           D4      D5
A2004 3235179 3985995

i.e.

> TRIANGLE2010["A2004","D5"]-TRIANGLE2010["A2004","D4"]
[1] 750816

Pour le montant payé au total en 2005, on utilise

> CUMULPAIEMENTS2010=TRIANGLE2010
> INCREMENTPAIEMENTS2010=CUMULPAIEMENTS2010
> INCREMENTPAIEMENTS2010[,2:9]=CUMULPAIEMENTS2010[,2:9]-
+ CUMULPAIEMENTS2010[,1:8]
> sum(diag(as.matrix(INCREMENTPAIEMENTS2010[1:4,4:1])))
[1] 2729241

On attaque maintenant la méthode Chain Ladder. Les coefficients de développement sont

> lambda=rep(NA,8)
> for(k in 1:8){
+ lambda[k]=(sum(TRIANGLE2010[1:(9-k),k+1])/
+ sum(TRIANGLE2010[1:(9-k),k]))}
> lambda[1:5]
[1] 3.474193 1.704149 1.460866 1.161765 1.095763

On peut alors obtenu ce qui sera payé en 2011 pour les sinistres survenus en 2010,

> TRIANGLE2010["A2010","D1"]*(lambda[1]-1)
[1] 931993.8

alors que pour le montant payé en 2011 pour les sinistres survenus en 2009,

> TRIANGLE2010["A2009","D2"]*(lambda[2]-1)
[1] 1000686

Pour le montant payé en 2010 et 2011 pour les sinistres survenus en 2010,

> TRIANGLE2010["A2010","D1"]*lambda[1]*lambda[2]
[1] 2230186

Pour les sinistres survenus en 2008, la charge total estimée sera

> TRIANGLE2010["A2008","D3"]*prod(lambda[3:8])
[1] 5531128

ce qui fait un montant de provision

> TRIANGLE2010["A2008","D3"]*(prod(lambda[3:8])-1)
[1] 3111267

On rajoute maintenant une diagonale dans notre triangle, i.e. une nouvelle année de paiements (en 2011).
Le montant payé en 2011 pour les sinistres survenus en 2010 est

> TRIANGLE2011["A2010","D2"]-TRIANGLE2011["A2010","D1"]
[1] 986608

à comparer avec

> TRIANGLE2010["A2010","D1"]*(lambda[1]-1)
[1] 931993.8

prédit par la méthode Chain-Ladder, soit une erreur de l’ordre de 5%,

> (TRIANGLE2011["A2010","D2"]-
+ TRIANGLE2011["A2010","D1"])/
+ (TRIANGLE2010["A2010","D1"]*(lambda[1]-1))-1
[1] 0.05859927

Si on compare maintenant les coefficients de transition de la méthode Chain Ladder,

> lambda2011=rep(NA,9)
> for(k in 2:9){
+ lambda2011[k]=(sum(TRIANGLE2011[1:(10-k),k+1])/
+ sum(TRIANGLE2011[1:(10-k),k]))}
> lambda2011[1:5]
[1]       NA 1.747333 1.457413 1.173852 1.103824
> lambda[1:5]
[1] 3.474193 1.704149 1.460866 1.161765 1.095763

Le coefficient permettant de passer de la colonne 2 à la colonne 3 augmente (légèrement). Si on regarde la prédiction du montant de paiements effectué en 2012 pour les sinistres survenus en 2009,

> TRIANGLE2011["A2009","D3"]*(lambda2011[3]-1)
[1] 1310258

qui sont à comparer avec les

> TRIANGLE2010["A2009","D2"]*lambda[2]*(lambda[3]-1)
[1] 1116132

estimés un an auparavant. La variation entre les estimations est de l’ordre de 17%,

> TRIANGLE2011["A2009","D3"]*(lambda2011[3]-1)/
+ (TRIANGLE2010["A2009","D2"]*lambda[2]*(lambda[3]-1))-1
[1] 0.1739278

Entre 2010 et 2012, on peut estimer que l’on payera

> TRIANGLE2011["A2010","D2"]*lambda2011[2]
[1] 2382128

pour les sinistres survenus en 2010.
Si on se focalise sur les sinistres survenus en 2008, la prédiction de la charge finale est ici

> TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9])
[1] 5660771

à comparer avec la prédiction faire en 2010,

> TRIANGLE2010["A2008","D3"]*prod(lambda[3:8])
[1] 5531128

i.e. notre nouvelle prédiction est 2% plus élevée,

> TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9])/
+ (TRIANGLE2010["A2008","D3"]*prod(lambda[3:8]))-1
[1] 0.02343867

Le provisions pour sinistres à payer pour les sinistres survenus en 2008 est estimé, fin 2011, à

> TRIANGLE2011["A2008","D4"]*(prod(lambda2011[4:9])-1)
[1] 2177641

alors que fin 2010, un montant de provision avait été constitué, à hauteur de

> TRIANGLE2010["A2008","D3"]*(prod(lambda[3:8])-1)
[1] 3111267

On notera que la différence

> TRIANGLE2011["A2008","D4"]*(prod(lambda2011[4:9])-1)-
+ TRIANGLE2010["A2008","D3"]*(prod(lambda[3:8])-1)
[1] -933626.7

vient des paiements effectués en 2011, et du changement dans l’estimation de charge ultime

> TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9])-
+ (TRIANGLE2010["A2008","D3"]*prod(lambda[3:8]))
[1] 129642.3
> TRIANGLE2011["A2008","D4"]-TRIANGLE2011["A2008","D3"]
[1] 1063269
> TRIANGLE2011["A2008","D4"]*prod(lambda2011[4:9])-
+ (TRIANGLE2010["A2008","D3"]*prod(lambda[3:8]))-
+ (TRIANGLE2011["A2008","D4"]-TRIANGLE2011["A2008","D3"])
[1] -933626.7

On attaque enfin la dernière partie, avec la régression de Poisson. La prédiction du montant de paiements faits en 2012 pour les sinistres survenus en 2009 est ici

> Y=as.vector(as.matrix(INCREMENTPAIEMENTS2011))
> ANNEE=rep(2002:2010,10)
> DEVEL=rep(1:10,each=9)
> baseTriangle=data.frame(Y,A=as.factor(ANNEE),
+ D=as.factor(DEVEL))
> reg=glm(Y~A+D,data=baseTriangle,famil=poisson(link="log"))
> nbase=data.frame(A="2009",D="4")
> predict(reg,newdata=nbase,type="response")
      1 
1310258

qui correspond à ce que donnait la méthode Chain Ladder

> TRIANGLE2011["A2009","D3"]*(lambda2011[3]-1)
[1] 1310258

(mais on savait que les deux modèles prédisaient la même chose). Pour les sinistres survenus en 2005, le montant de provisions à constituer est

> nbase=data.frame(A="2005",D=c("8","9","10"))
> predict(reg,newdata=nbase,type="response")
       1        2        3 
247190.0 370179.3  92268.5 
> sum(predict(reg,newdata=nbase,type="response"))
[1] 709637.8

qui, là encore, correspond à ce que prédisait la méthode Chain Ladder.

> TRIANGLE2011["A2005","D7"]*(prod(lambda2011[7:9])-1)
[1] 709637.8

Incréments négatifs dans les triangles de paiements

Considérons le triangle d’incréments de payements suivants,

> PAID
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 4372 4411 4428 4435 4456
[2,] 3367 4659 4696 4720 4730   NA
[3,] 3871 5345 5338 5420   NA   NA
[4,] 4239 5917 6020   NA   NA   NA
[5,] 4929 6794   NA   NA   NA   NA
[6,] 5217   NA   NA   NA   NA   NA
> INC=PAID
> INC[,2:6]=PAID[,2:6]-PAID[,1:5]
> INC
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,] 3209 1163   39   17    7   21
[2,] 3367 1292   37   24   10   NA
[3,] 3871 1474   -7   82   NA   NA
[4,] 4239 1678  103   NA   NA   NA
[5,] 4929 1865   NA   NA   NA   NA
[6,] 5217   NA   NA   NA   NA   NA

Comme on peut le voir, sur la troisième ligne, on a un incrément de paiementnégatif. A priori, ce n’est pas forcément gênant. En particulier, on peut faire tourner la méthode Chain Ladder,

> lambda=rep(NA,5)
> for(k in 1:5){
+ lambda[k]=(sum(PAID[1:(6-k),k+1])/
+             sum(PAID[1:(6-k),k]))}
> lambda
[1] 1.380933 1.008476 1.008515 1.001858 1.004735
> PROJECTION=PAID
> for(k in 1:5){
+ PROJECTION[((7-k):6),k+1]=
+ PROJECTION[((7-k):6),k]*lambda[k]}
> sum(PROJECTION[,6]-
+ diag(PROJECTION[,6:1]))
[1] 2469.703

et la sortie coïncide avec la fonction R,

> MackChainLadder(PAID)
MackChainLadder(Triangle = PAID)

Latest Dev.To.Date Ultimate    IBNR Mack.S.E CV(IBNR)
1  4,456       1.000    4,456     0.0    0.000      NaN
2  4,730       0.995    4,752    22.4    0.146  0.00652
3  5,420       0.993    5,456    35.8    2.405  0.06721
4  6,020       0.985    6,111    91.3   41.679  0.45629
5  6,794       0.977    6,956   161.5   71.620  0.44334
6  5,217       0.707    7,376 2,158.6   95.750  0.04436

Totals
Latest:            32,637.00
Dev:                    0.93
Ultimate:          35,106.70
IBNR:               2,469.70
Mack S.E.:            146.62
CV(IBNR):  0.059366227164502
Message d'avis :
In Mack.S.E(CL[["Models"]], FullTriangle, est.sigma = est.sigma,
'loglinear' model to estimate sigma_n doesn't appear appropriate
p-value > 5.
est.sigma will be overwritten to 'Mack'.
Mack's estimation method will be used instead.

On notera le message d’avis, laissant entendre qu’il peut y avoir un petit soucis. En fait, le soucis apparaît clairement si on souhaite faire une régression de Poisson. Car autant les valeurs non-entières ne sont pas trop gênantes, autant les valeurs négatives le sont !

> Y=as.vector(INC)
> D=rep(1:6,each=6)
> A=rep(2001:2006,6)
> base=data.frame(Y,D,A)
> reg=glm(Y~as.factor(D)+as.factor(A),
+     data=base,family=poisson(link="log"))
Erreur dans eval(expr, envir, enclos) :
les valeurs négatives sont interdites pour la famille poisson

Il faut alors trouver une solution si on a des incréments négatifs. Et aucune solution ne sera intellectuellement satisfaisante, car avec une valeur négative, il n’y a aucune légitimité pour continuer à utiliser un modèle Poissonnien. Donc on va bricoler…

  • Prendre de l’argent à droite et à gauche

La première solution consiste à dire que cet incrément n’a pas de raison d’être. On va donc aller chercher de l’argent dans la colonne avant (sur la même ligne, on va supposer qu’il s’agit d’un problème sur les cadences de paiements) ou sur la colonne après. Au lieu d’avoir

[3,] 3871 1474   -7   82   NA   NA

On peut renflouer en prenant à gauche,

[3,] 3871 1467   0    82   NA   NA

ou à droite

[3,] 3871 1474   0    75   NA   NA

En faisant ces opérations, on change seulement localement la courbe de paiements pour la troisième année de survenance. On peut se demander si le fait de prendre à gauche a un impact sur l’estimation du montant de provisions, ou sur l’incertitude associée. Pour ça, on peut utiliser les fonctions suivantes,

library(ChainLadder)
CL=function(T){
sum(
MackChainLadder(T)$FullTriangle[,ncol(T)]
-diag(T[,ncol(T):1]))}
CL.SE=function(T){
MackChainLadder(T)$Total.Mack.S.E}
CumInc=function(T){
m=T
for(i in 2:nrow(T)){m[,i]=apply(T[,1:i],1,sum)}
return(m)}

On peut alors regarder, si on prend à droite ou à gauche ce qui se passe, ou ce qui se passe si on renfloue non pas à 0 mais à une valeur strictement positive (e.g. 1),

VCL=rep(NA,8)
for(j in 1:8){
INCd=INC
INCd[3,2:4]=INC[3,2:4]+c(1-j,+7,j-1-7)
VCL[j]=CL(CumInc(INCd))}
plot(1:8,VCL,type="b",col="blue",xlim=c(0,9),ylim=c(2462,2474))
VCL=rep(NA,10)
for(j in 1:10){
INCd=INC
INCd[3,2:4]=INC[3,2:4]+c(1-j,+9,j-1-9)
VCL[j]=CL(CumInc(INCd))}
lines(0:9,VCL,type="b",pch=0,col="red")

avec la courbe bleu si on renfloue à 0, et rouge si on renfloue à 1, pour l’estimation du montant total de provisions,

VCL=rep(NA,8)
for(j in 1:8){
INCd=INC
INCd[3,2:4]=INC[3,2:4]+c(1-j,+7,j-1-7)
VCL[j]=CL.SE(CumInc(INCd))}
plot(1:8,VCL,type="b",col="blue",xlim=c(0,9),ylim=c(130,145))
VCL=rep(NA,10)
for(j in 1:10){
INCd=INC
INCd[3,2:4]=INC[3,2:4]+c(1-j,+9,j-1-9)
VCL[j]=CL.SE(CumInc(INCd))}
lines(0:9,VCL,type="b",pch=0,col="red")

avec cette fois l’impact des transferts sur l’écart-type. En abscisse, on a (à peu de choses près) le montant enlevé sur la colonne de droite. Bref, les transferts venant de la gauche ou de la droite sont une solution, mais le choix d’où vient l’argent ne sera pas neutre, ni sur l’estimation, ni sur la variance de l’estimation (et l’intervalle de confiance).

  • Jouer à faire des translations…

Une autre piste pourrait être de noter que, dans le modèle linéaire, si on translate nos données (vers le haut), on ne change pas la pente, et la constante est juste augmenté (de la taille de la translation).

> lm(dist~speed,data=cars)

Call:
lm(formula = dist ~ speed, data = cars)

Coefficients:
(Intercept)        speed
-17.579        3.932

> lm((dist+10)~speed,data=cars)

Call:
lm(formula = (dist + 10) ~ speed, data = cars)

Coefficients:
(Intercept)        speed
-7.579        3.932

Autrement dit, la prédiction faite par notre modèle n’est pas modifiée par la translation, à condition d’opérer la translation ensuite sur la prédiction. Sur le dessin ci-dessous, on fait pareil, mais avec une régression de Poisson: la courbenoire est sur les données brutes, la courbe bleu est obtenue sur les points translatés vers le haut, et la rouge si on translate la prédiction sur les points translatés,
http://freakonometrics.blog.free.fr/public/perso4/glm-translation.gif
Si la translation est trop importante, la prédiction finale (en rouge) se retrouve davantage éloignée de la vraie valeur (en noir).
Pareillement, on peut translater nos paiements vers le haut dans les triangles, de manière à avoir des paiements positifs. On peut commencer par translater tout le triangle, par exemple en translatant juste assez pour que tous les incréments soient positifs ou nuls,

> k=7
> Y=as.vector(INC)+k
> D=rep(1:6,each=6)
> A=rep(2001:2006,6)
> base=data.frame(Y,D,A)
> reg=glm(Y~as.factor(D)+as.factor(A),
+ data=base,family=poisson(link="log"))
> Yp=predict(reg,type="response",
+ newdata=base)-k
> sum(Yp[is.na(Y)==TRUE])
[1] 2508.620

On peut aussi se demander ce qui se passerait si on décalait davantage vers le haut. On l’a vu sur l’animation, plus la translation est importante, plus la prédiction se décale. Un stratégie pourrait être de faire une estimation pour plusieurs valeurs de translations – de telle sorte que tous les incréments soient positifs ou nuls – et ensuite d’extrapoler pour savoir ce qui se passerait si on ne translatait pas,

> K=7:20
> R=rep(NA,length(K))
> for(i in 1:length(K)){
+ k=K[i]
+ Y=as.vector(INC)+k
+ D=rep(1:6,each=6)
+ A=rep(2001:2006,6)
+ base=data.frame(Y,D,A)
+ reg=glm(Y~as.factor(D)+as.factor(A),
+ data=base,family=poisson(link="log"))
+ Yp=predict(reg,type="response",
+ newdata=base)-k
+ R[i]=sum(Yp[is.na(Y)==TRUE])}
> plot(K,R,xlim=c(0,20),ylim=c(2465,max(R)))
> abline(lm(R~K),col="blue")
> (yp=predict(lm(R~K),newdata=(K=0)))
1
2470.199
> points(0,yp,col="red",pch=19)

 

On a cette fois une estimation plus proche de celle que nous avions en utilisant directement, sur le triangle brut, la méthode chain ladder. Mais on peut aller un peu plus loin. Car dans le triangle, on a translaté toutes les valeurs, alors que seule une était négative. On pourrait par exemple translater seulement les valeurs de la troisième colonne,

> K=7:20
> R=rep(NA,length(K))
> for(i in 1:length(K)){
+ k=K[i]
+ Y=as.vector(INC)
+ D=rep(1:6,each=6)
+ A=rep(2001:2006,6)
+ Y[D==3]=Y[D==3]+k
+ base=data.frame(Y,D,A)
+ reg=glm(Y~as.factor(D)+as.factor(A),
+ data=base,family=poisson(link="log"))
+ Yp=predict(reg,type="response",
+ newdata=base)
+ Yp[D==3]=Yp[D==3]-k
+ R[i]=sum(Yp[is.na(Y)==TRUE])}
> predict(lm(R~K),newdata=(K=0))
1
2469.703

Tiens, on retombe exactement sur l’estimateur de la méthode chain ladder… Étonnant, non ?

Une petite précision. Dans la vraie vie, on ne devrait pas avoir d’incréments négatifs dans les triangles de paiements. Maintenant, il faut reconnaître qu’on en voit apparaître lorsque l’on bootstrappe les résidus, et que l’on génère de pseudo-triangles…
Si des praticiens ont des commentaires sur les incréments négatifs, les commentaires sont ouverts (et peu modérés), donc racontez nous comment vous faites, je suis preneur…