Chain Ladder, avec R

Un billet rapide pour mettre en ligne des parties du code tapé en cours, mercredi dernier. On avait commencé par convertir la feuille du classeur excel en un fichier texte, pour faciliter la lecture,

> setwd("C:\\Users\\savsalledecours\\Desktop")
> triangle=read.table("exACT2040.csv",header=TRUE,sep=";")	
> triangle
  ANNEE   X0   X1   X2   X3   X4   X5
1  2000 3209 4372 4411 4428 4435 4456
2  2001 3367 4659 4696 4720 4730   NA
3  2002 3871 5345 5398 5420   NA   NA
4  2003 4239 5917 6020   NA   NA   NA
5  2004 4929 6794   NA   NA   NA   NA
6  2005 5217   NA   NA   NA   NA   NA

L’idée – quand on importe un triangle – est de récupérer une base sous la forme précédente, avec des valeurs manquantes dans la partie inférieure du triangle (on verra l’intérêt quand on fait une régression). On avait ensuite calculé les facteurs de transition, et en même temps complété le triangle,

> T=triangle[,2:7]
> rownames(T)=triangle$ANNEE
> T2=T
> n=ncol(T)
> L=rep(NA,n-1)
> for(j in 1:(n-1)){
+ L[j]=sum(T[1:(n-j),j+1])/sum(T[1:(n-j),j])
+ T2[(n-j+1):n,j+1]=L[j]*T2[(n-j+1):n,j]
+ }

Les facteurs de transition sont ici,

> L
[1] 1.380933 1.011433 1.004343 1.001858 1.004735

et le triangle complété

> T2
       X0       X1       X2       X3       X4       X5
2000 3209 4372.000 4411.000 4428.000 4435.000 4456.000
2001 3367 4659.000 4696.000 4720.000 4730.000 4752.397
2002 3871 5345.000 5398.000 5420.000 5430.072 5455.784
2003 4239 5917.000 6020.000 6046.147 6057.383 6086.065
2004 4929 6794.000 6871.672 6901.518 6914.344 6947.084
2005 5217 7204.327 7286.691 7318.339 7331.939 7366.656

Le montant de provision est ici en faisant la différence entre la charge ultime (dans la dernière colonne) et les derniers paiements observés (sur la seconde diagonale)

> CU=T2[,n]
> Pat=diag(as.matrix(T2[,n:1]))
> Ri=CU-Pat
> R=sum(Ri)

soit, numériquement

> R
[1] 2426.985

On avait alors vu que l’on pouvait calculer un tail factor, en supposant une décroissance exponentielle des facteurs de transition, et on rajoutait alors une colonne correspondant au montant ultime, par année d’accident,

> logL=log(L-1)
> t=1:5
> b=data.frame(logL,t)
> reg=lm(logL~t,data=b)
> logLp=predict(reg,newdata=data.frame(t=6:100))
> Lp=exp(logLp)+1
> Linf=prod(Lp)
> T3=T2
> T3$Xinf=T3$X5*Linf

On a ici

> T3
       X0       X1       X2       X3       X4       X5     Xinf
2000 3209 4372.000 4411.000 4428.000 4435.000 4456.000 4459.149
2001 3367 4659.000 4696.000 4720.000 4730.000 4752.397 4755.755
2002 3871 5345.000 5398.000 5420.000 5430.072 5455.784 5459.639
2003 4239 5917.000 6020.000 6046.147 6057.383 6086.065 6090.366
2004 4929 6794.000 6871.672 6901.518 6914.344 6947.084 6951.993
2005 5217 7204.327 7286.691 7318.339 7331.939 7366.656 7371.862

(je laisse reprendre le code pour calculer le montant de provisions). Enfin, on avait montré comment utiliser une régression pondérée, pour calculer les facteurs de transition,

> T4=as.matrix(T$X0,n,1)
> for(j in 1:(n-1)){
+ Y=T[,j+1]
+ X=T[,j]
+ base=data.frame(X,Y)
+ reg=lm(Y~0+X,weights=1/X)
+ T4=cbind(T4,
+ predict(reg,
+ newdata=data.frame(X=T4[,j]
+ )))
+ }

Ce qui donnait la même projection que la méthode Chain Ladder

> T4
  [,1]     [,2]     [,3]     [,4]     [,5]     [,6]
1 3209 4431.414 4482.076 4501.543 4509.909 4531.263
2 3367 4649.601 4702.758 4723.184 4731.961 4754.367
3 3871 5345.591 5406.705 5430.188 5440.279 5466.039
4 4239 5853.775 5920.698 5946.414 5957.464 5985.673
5 4929 6806.619 6884.435 6914.337 6927.186 6959.986
6 5217 7204.327 7286.691 7318.339 7331.939 7366.656

La suite mercredi prochain, même si on risque d’aller très vite sur la méthode de Mack (et les calculs d’erreur quadratique moyenne pour arriver à la régression de Poisson). A suivre donc…

Easter

This morning, there was an interesting post entitled “why does Easter move around so much?” online on http://economist.com/blogs/economist-explains/…

In my time series classes, I keep saying that sometimes, series can exhibit seasonlity, but the seasonal effect can be quite irregular. It is the cas for river levels, where snowmelt can have a huge impact, and it is irregular. Similarly, chocolate sales (even monthly, or quarterly) depends on Easter. Because it can be either in March, or in April, the seasonal pattern is not as regular as flower sales for instance (Valentine beeing always on February 14th, as far as I remember). If we look at the word eggs on http://google.com/trends/q=eggs…, we do observe a cycle related to Easter.

The title of the article published by http://economist.com/blogs/economist-explains/… claims that there is a lot of variability on Eater’s day. Let us check ! The answer to the question “When is Easter ?” can be the following (if we want a short answer): Easter Sunday is the first Sunday after the first full moon after vernal equinox. For more details, see e.g. http://ortelius.de/east. The algorithm used to compute the date of Easter can is online, on http://smart.net/~mmontes/….

> century = year/100
> G = year % 19
> K = (century - 17)/25
> I = (century - century/4 - (century - K)/3 + 19*G + 15) % 30
> I = I - (I/28)*(1 - (I/28)*(29/(I + 1))*((21 - G)/11))
> J = (year + year/4 + I + 2 - century + century/4) % 7
> L = I - J
> EasterMonth = 3 + (L + 40)/44
> EasterDay = L + 28 - 31*(EasterMonth/4)

Actually, this algorithm can be found in some R packages. Here we use the date of Easter from AD 1000 and AD 3000,

> library(timeDate)
> E=Easter(1000:3000)
> D=as.Date(E)
> table(months(D))/2001

    april     march 
0.7651174 0.2348826

(April being before March, in the alphabetical order) If we look at the distribution of the date, it is the following, the starting point being March 1st,

> J=as.numeric(D-as.Date(paste("01/03/",1000:3000,sep=""),"%d/%m/%Y"))
> hist(J,breaks=seq(20,55),col="light green")

And if we look at the autocorrelation function, we can observe that indeed, after 19 years, there is a strong correlation (that could be seen in the algorithm given previously),

> plot(acf(J))

But in order to get a better understanding of the dynamics, we can also look at transiftion matrices. Define

> Q=quantile(J,seq(0,1,by=.25))
> Q[1]=Q[1]-1
> C=cut(J,Q)

Then, the one year transition matrix is (in %)

> k=1; n=length(C)
> B=data.frame(X1=(C[1:(n-k)]),X2=(C[(k+1):n]))
> (T=table(B$X1,B$X2))

          (20,31] (31,39] (39,46] (46,55]
  (20,31]       0       0     265     277
  (31,39]     316       0      13     182
  (39,46]     224     264       0       0
  (46,55]       1     247     211       0
> P=T/apply(T,1,sum)
> round(P*1000)/10

          (20,31] (31,39] (39,46] (46,55]
  (20,31]     0.0     0.0    48.9    51.1
  (31,39]    61.8     0.0     2.5    35.6
  (39,46]    45.9    54.1     0.0     0.0
  (46,55]     0.2    53.8    46.0     0.0

I.e. if  Easter was early in the year (say in March, in the first quartile), then very likeliy, the year after, it will be late in the year (with 50% chance in the third quartile, and 50% chance in the fourth one).

Le libre, et l’universitaire

On a pu lire ces derniers jours plusieurs articles sur le libre accès aux publications scientifique. Je pense à la tribune “qui a peur de l’open access publiée sur http://lemonde.fr/sciences/… . Etant hébergé par http://hypotheses.org/ qui milite en faveur du libre accès, et d’une plus vaste diffusion des travaux et de la culture scientifique, je n’y suis pas insensible. Sur les publications, la conférence des présidents d’université a aussi publié un document intéressant, intitulé “bien universel, par essence, l’article scientifique, n’est pas un bien comme les autres” en ligne sur http://cpu.fr/

En fait, le but de mon billet aujourd’hui était non pas de discuter sur le libre accès aux publications scientifiques, mais à revenir sur l’utilisation des logiciels libres dans le monde universitaire, et plus spécifiquement sur la place du logiciel libre dans l’enseignement.

En lisant Scott Wilson, qui a publié “open source in higher education: how far have we come?”,  sur http://guardian.co.uk/higher-education-network/… je l’ai trouvé un peu trop optimiste. Dans son article  il note que “universities are ahead of the curve in adopting open source, we should now lead the public sector in exploring its full potential“. Jusqu’à aujourd’hui (ou disons récemment pour être honnête), j’étais assez étranger au débat. Certes, j’utilise des logiciels libres, mais juste parce que je les trouve meilleurs. Ce n’était pas du militantisme, juste du pragmatisme. Et depuis quelque temps, je prends conscience que j’ai (ou que nous avons, en tant qu’universitaires) une responsabilité collective. Par exemple dans mes enseignements, j’utilise R (autant le nommer, ce n’est un secret pour personne que la moitié des billets sur mon blog contiennent du R). Mais quand j’étais étudiant, j’ai appris à manipuler des logiciels de statistiques payant (appelons les S pour utiliser une notation générique qui devrait m’éviter des procès pour diffamation). J’ai aussi utilisé S lorsque j’ai travaillé en tant qu’actuaire, car c’était installé sur nos machines, et parce que j’avais appris ce logiciel. Sans parler de premières amours, disons qu’on garde une forme d’attachement, ou de sécurité, à travailler avec les logiciels que l’on a vus (voire appris) durant nos études. On a peur de changer de logiciel, on aime la stabilité
(car admettons le, on est tous aussi nuls les uns que les autres avec l’informatique). C’est un peu ce que l’on peut lire dans un article paru il y a quelques jours dans http://ledevoir.com/politique/quebec/… où l’on apprend que “Québec adopte un autre décret pour la mise à jour de 76 000 postes de la fonction publique, sans appel d’offres, par des logiciels vendus par […]. La facture totale est de 30 millions, pour les licences seulement. Cette mesure d’exception vise à « assurer la sécurité de nos postes »” (c’est moi qui met en gras). Ce n’est pas de sécurité informatique face à du piratage dont on parle ici, mais de sécurité du personnel qui utilise lesdits logiciels (le confort en quelque sorte). On comprends ainsi que S (et je mets dans le même panier tous les logiciels commerciaux, ou presque) offre des licences gratuites aux étudiants. On retrouve là la logique des fabricants de cigarettes il y a quelques décennies (suis-je optimiste ?) lorsqu’ils essayaient d’attirer les adolescents vers leurs cigarettes, tenant compte ainsi de l’addiction qui se crée ?

Comme je le disais au début de mon paragraphe précédent c’est par pragmatisme que, pour ma recherche, je me suis tourné vers les logiciels libres. Mais pour mon enseignement, je n’avais jamais vraiment réfléchi de manière profonde (parce que je ne fais pas de cours d’informatique: si on utilise un langage c’est parce que l’outil informatique est indispensable pour faire de la modélisation). Par simplicité (ou par fainéantise), j’enseignais des langages que je connaissais. J’utilise R dans mes cours, mais parmi mes proches voisins de bureau, Jean-Philippe utilise S (pour faire des régressions avancées en assurance non-vie) et Mathieu utilise un tableur (que l’on pourrait appeler T, pour…. faire comprendre qu’un tableur suffit pour faire de l’actuariat ? que l’actuariat est plus proche de la comptabilité que de la modélisation statistique ?). Récemment, plusieurs petites phrases m’ont poussé à m’interroger. Je pense que ça a commencé quand un(e) étudiant(e) a écrit dans l’évaluation du cours de modèle de prévision que je donnais voilà quelque temps, que mon cours se résumait à”faire de la saisi de données dans un programme gratuit sur internet” (je passe le fait que dans les questions posées aux étudiants qui remplissent l’évaluation, il est demandé ce que les étudiants pensent des présentations faites avec le logiciel P, comme si faire des présentation en LaTeX n’était pas admissible). Il y a eu ce commentaire l’autre jour sur mon blog, où quelqu’un me demandait le code en S pour refaire ce que je faisais en R. Il y a eu ce questionnement l’autre jour : “est-ce gênant si les étudiants du baccalauréat en actuariat sortent de la formation sans connaître S ?” (que j’avais entendu à Rennes 1 lors des discussions du contenu des cours du Master de statistique et d’économie “il est bon que nos étudiants aient 20 heures de cours à S, ça leur permet de mieux se placer sur le marché du travail“).

Le modèle économique est assez incroyable quand on y pense : avec des fonds publics, l’université offre un service avant-vente à des entreprises qui vendent des logiciels commerciaux. Alors qu’il existe du libre qui fait mieux. Et qui continuera à faire mieux si la communauté d’utilisateurs est active. Je pense que les logiciels commerciaux n’ont pas leur place dans les formations universitaires, et qu’il serait temps que les universitaires prennent leur part de responsabilité dans ce qui s’est mis en place Je ne souhaite pas que les universitaires prennent leur part du profit même si ça serait la moindre des choses (il suffit de survoler la page comparison of statistical packages sur wikipedia pour apprendre que rendre accro un étudiant à un logiciel qu’il utilisera pendant 10 ans, c’est assurer entre 10,000 $ et 60,000 $ de licence (pour un utilisateur), sans parler des formations, des livres, etc). Je pense qu’il y va de notre responsabilité en tant qu’universitaire de promouvoir le logiciel libre. Cela permet déjà de comprendre ce qu’on utilise (je faisais un billet la semaine passée, suite à des soucis rencontrés sur le tableur – commercial – T), mais en plus, comme le notait Scott Wilson, “open source not only promotes creativity but helps make technology more democratic allowing a community to work together to solve common problems.” Et je pense qu’il a raison. Maintenant, je suis nul pour faire des tribunes, et je voulais juste que mes collègues prennent conscience de certains aspects d’un débat qui ne fera que gonfler dans les mois qui viennent…

Benford law and lognormal distributions

Benford’s law is nowadays extremely popular (see e.g. http://en.wikipedia.org/…). It is usually claimed that, for a given set data set, changing units does not affect the distribution of the first digit. Thus, it should be related to scale invariant distributions. Heuristically, scale (or unit) invariance means that the density of the measure https://latex.codecogs.com/gif.latex?%20X (or probability function) https://latex.codecogs.com/gif.latex?f(x) should be proportional to https://latex.codecogs.com/gif.latex?f(kx). Thus, because densities integrate to 1, the proportionality coefficient has to be https://latex.codecogs.com/gif.latex?k^{-1}, and therefore, https://latex.codecogs.com/gif.latex?f should satisfy the following functional equation, https://latex.codecogs.com/gif.latex?%20kf(kx)=f(x), for all https://latex.codecogs.com/gif.latex?%20x in https://latex.codecogs.com/gif.latex?%20(1,\infty) and https://latex.codecogs.com/gif.latex?%20k in https://latex.codecogs.com/gif.latex?%20(0,\infty). The solution of this functional equation is https://latex.codecogs.com/gif.latex?%20f(x)=x^{-1}, I guess this can be proved easily solving ordinary differential equation

https://latex.codecogs.com/gif.latex?%20\frac{d}{dk}%20(kf(kx))=0

Now if https://latex.codecogs.com/gif.latex?%20D denotes the first digit of https://latex.codecogs.com/gif.latex?%20X, in base 10, then

https://latex.codecogs.com/gif.latex?%20\mathbb{P}(D=d)=\frac{\displaystyle{\int_d^{d+1}%20f(x)dx}}{{\displaystyle{\int_1^{10}%20f(x)dx}}}=\cdots=\frac{\displaystyle{\log\left(1+\frac{1}{d}\right)}}{\log(10)}Which is the so-called Benford’s law. So, this distribution looks like that

> (benford=log(1+1/(1:9))/log(10))
[1] 0.30103000 0.17609126 0.12493874 0.09691001 0.07918125 
[6] 0.06694679 0.05799195 0.05115252 0.04575749
> names(benford)=1:9
> sum(benford)
[1] 1
> barplot(benford,col="white",ylim=c(-.045,.3))
> abline(h=0)

To compute the empirical distribution from a sample, use the following function

> firstdigit=function(x){
+ if(x>=1){x=as.numeric(substr(as.character(x),1,1)); zero=FALSE}
+ if(x<1){zero=TRUE}
+ while(zero==TRUE){
+ x=x*10; zero=FALSE
+ if(trunc(x)==0){zero=TRUE}
+ }
+ return(trunc(x))
+ }

and then

> Xd=sapply(X,firstdigit)
> table(Xd)/1000

In Benford’s Law: An Empirical Investigation and a Novel Explanation, we can read

It is not a mathematical article, so do not expect any formal proof in this paper. At least, we can run monte carlo simulation, and see what’s going on if we generate samples from a lognormal distribution with variance https://latex.codecogs.com/gif.latex?%20\sigma^2. For instance, with a unit variance,

> set.seed(1)
> s=1
> X=rlnorm(n=1000,0,s)
> Xd=sapply(X,firstdigit)
> table(Xd)/1000
Xd
    1     2     3     4     5     6     7     8     9 
0.288 0.172 0.121 0.086 0.075 0.072 0.073 0.053 0.060 
> T=rbind(benford,-table(Xd)/1000)
> barplot(T,col=c("red","white"),ylim=c(-.045,.3))
> abline(h=0)

Clearly, it not far away from Benford’s law. Perhaps a more formal test can be considered, for instance Pearson’s https://latex.codecogs.com/gif.latex?%20\chi^2 (goodness of fit) test.

> chisq.test(T,p=benford)

	Chi-squared test for given probabilities

data:  T 
X-squared = 10.9976, df = 8, p-value = 0.2018

So yes, Benford’s law is admissible ! Now, if we consider the case where https://latex.codecogs.com/gif.latex?%20\sigma is smaller (say 0.9), it is a rather different story,

compared with the case where https://latex.codecogs.com/gif.latex?%20\sigma is larger (say 1.1)

It is possible to generate several samples (always the same size, here 1,000 observations), just change the variance parameter https://latex.codecogs.com/gif.latex?%20\sigma and compute the https://latex.codecogs.com/gif.latex?%20p-value of the test. There might be one tricky part: when generating samples from lognormal distributions with small variance, it might be possible that some digits do not appear at all. On that case, there is a problem with the test. So we just use here

> T=table(Xd)
> T=T[as.character(1:9)]
> T[is.na(T)]=0
> PVAL[i]=chisq.test(T,p=benford)$p.value

Boxplots of the https://latex.codecogs.com/gif.latex?%20p-value of the test are the following,

When https://latex.codecogs.com/gif.latex?%20\sigma is too small, it is clearly not Benford’s distribution: for half (or more) of our samples, the https://latex.codecogs.com/gif.latex?%20p-value is lower than 5%. On the other hand, when https://latex.codecogs.com/gif.latex?%20\sigma is large (enough), Benford’s distribution is the distribution of the first digit of lognormal samples, since 95% of our samples have  https://latex.codecogs.com/gif.latex?%20p-values higher than 5% (and the distribution of the https://latex.codecogs.com/gif.latex?%20p-value is almost uniform on the unit interval). Here is the proportion of samples where the https://latex.codecogs.com/gif.latex?%20p-value was lower than 5% (on 5,000 generations each time)

Note that it is also possible to compute the https://latex.codecogs.com/gif.latex?%20p-value of Komogorov-Smirnov test, testing if the https://latex.codecogs.com/gif.latex?%20p-value has a uniform distribution,

> ks.test(PVAL[,s], "punif")$p.value

Indeed, if https://latex.codecogs.com/gif.latex?%20\sigma is larger than 1.15 (around that value), it looks like Benford’s law is a suitable distribution for the first digit.

Rationality, and MS Excel (and other calculators)

This morning, Mathieu had a nice experience in his course on computational method in actuarial science. But let us start with some mathematical formal definitions.

First, recall that https://latex.codecogs.com/gif.latex?y^x is – somehow – a standard expression. No one should be surprised to see such an expression. Generally (as explained in http://en.wikipedia.org/… ), this function is defined only when https://latex.codecogs.com/gif.latex?y\in\mathbb{R}_+. The idea is that the definition of https://latex.codecogs.com/gif.latex?y^x is that

https://latex.codecogs.com/gif.latex?y^x%20=%20\exp\left(x\log[y]\right)

And it is a definition. Such a function exists only if https://latex.codecogs.com/gif.latex?y\in\mathbb{R}_+ (maybe excluding https://latex.codecogs.com/gif.latex?0). This would be a standard definition in real-analysis.

Now, this ‘power’ function appears also in complex analysis, when dealing with unit roots. From instance, if  https://latex.codecogs.com/gif.latex?z=y^{\frac{1}{k}}e^{i%20\frac{2n\pi}{k}}, where https://latex.codecogs.com/gif.latex?y\in\mathbb{R}_+ and https://latex.codecogs.com/gif.latex?k\in\mathbb{N}_\star, for some https://latex.codecogs.com/gif.latex?n\in\mathbb{N}, then https://latex.codecogs.com/gif.latex?z^k=y. Thus, in complex-analysis it might be more complex to define properly https://latex.codecogs.com/gif.latex?y^x since it might not be unique. But we can relate (sometimes, when https://latex.codecogs.com/gif.latex?x is the inverse of an integer, or maybe a rational number ?) with roots of polynomial functions. So far, nothing new…

Let us get back to Mathieu’s problem. Actually, in his course, he wanted to compute https://latex.codecogs.com/gif.latex?(-8)^{\frac{1}{3}}. With a French version of Excel, entering

you do get https://latex.codecogs.com/gif.latex?-2. If you look at the ‘help’ window, you have some more details

It looks like this hat function can be used to define objects such as https://latex.codecogs.com/gif.latex?y^x. But with

you get

(meaning that this is a problem…). It is also possible to use the power (puissance in French) function of Excel,

Here, you also get

The weird part here is that, in the ‘help’ window, you can read that this power function can be used with any number in https://latex.codecogs.com/gif.latex?\mathbb{R}.

Another point… what about  ? Somehow, it is just the square of the previous one (with the fraction)… Here, typing

you get

(similarly with the power function). So clearly, it is not that simple to use this power function. Now, if you use Google (which is now my new online calculator when I am in class, when I cannot use R), if the power is a fraction (or to be more specific the inverse of an integer), then it works as Excel

 

you get

 But if you type (which should be close, from a continuity property of the power function)

 

you get

and similarly

On Wolfram Mathworld, enter

Mathematica does recognize that we try to deal with unit roots: the result is here

with – as expected – a numerical approximation

With Matlab, Mathieu did obtain the same as Mathematica (its decimal approximation). And to conclude, with R, Mathieu did obtain

> (-8)^(1/3)
[1] NaN
> (-8)^(.333333333333333)
[1] NaN

So for R, you cannot use this hat function on negative numbers.

Now, how can we interpret those outputs ?

1) My understanding is that clearly, with MS Excel, https://latex.codecogs.com/gif.latex?x^{ab}\neq%20\left(x^a\right)^bsince

https://latex.codecogs.com/gif.latex?(-8)^{\frac{2}{3}}\neq%20\left((-8)^{\frac{1}{3}}\right)^2

which is problematic. For instance, in insurance, with monthly discounts, we do have functions like https://latex.codecogs.com/gif.latex?u^{\frac{k}{12}}. What if

https://latex.codecogs.com/gif.latex?u^{\frac{k}{12}}\neq%20\left(u^{\frac{1}{12}}\right)^k

2) The problem comes – probably (MS Excel is not an open software, so it might be hard to check) –  from the fact that https://latex.codecogs.com/gif.latex?y^{\frac{1}{n}} is interpreted as an inverse of a (possibly) bijective function. To be more specific, https://latex.codecogs.com/gif.latex?x=y^{\frac{1}{n}} means that https://latex.codecogs.com/gif.latex?x^n=y. When https://latex.codecogs.com/gif.latex?n is an odd integer, then (in real-analysis) there is a unique inverse, and thus, https://latex.codecogs.com/gif.latex?y^{\frac{1}{n}} is uniquely defined, since https://latex.codecogs.com/gif.latex?x\mapsto%20x^n is a bijective https://latex.codecogs.com/gif.latex?\mathbb{R}\rightarrow\mathbb{R} function. This is what MS Excel (and Google) is doing: https://latex.codecogs.com/gif.latex?x\mapsto%20x^3 is a bijective https://latex.codecogs.com/gif.latex?\mathbb{R}\rightarrow\mathbb{R} function, so https://latex.codecogs.com/gif.latex?(-8)^{\frac{1}{3}} means that we need to find the unique (real) value https://latex.codecogs.com/gif.latex?x such that https://latex.codecogs.com/gif.latex?x^3=-8. Thus, somehow, it makes sense to return https://latex.codecogs.com/gif.latex?-2.

3) There is still a problem with Google, and Mathematica. That is fine to return unit roots in https://latex.codecogs.com/gif.latex?\mathbb{C}. But how comes there is only one value ? I mean, yes https://latex.codecogs.com/gif.latex?1+\sqrt{3}%20\%20i is a possible answer, since

https://latex.codecogs.com/gif.latex?(1+\sqrt{3}%20\%20i)^3=-8

but one can also observe that , and similarly, https://latex.codecogs.com/gif.latex?(-2)^3=-8 and

https://latex.codecogs.com/gif.latex?(1-\sqrt{3}%20\%20i)^3=-8

One can check with

With R, since we do not deal with power function here, but with roots, if we want to find https://latex.codecogs.com/gif.latex?x such that https://latex.codecogs.com/gif.latex?x^3=-8, the function is

> polyroot(c(8,0,0,1))
[1]  1+1.732051i -2+0.000000i  1-1.732051i

Which is different… Weird isn’t it ?

Solvabilité et provisionnement

Mercredi, nous allons aborder en cours les aspects de solvabilité des compagnies d’assurance IARD. Plus particulièrement, nous parlerons des provisions pour sinistres à payer, ou “provision for claims outstanding (PCO)” selon la terminologie anglaise, i.e. “the estimated total cost of ultimate settlement of all claims incurred before the date of record, whether reported or not, less any amounts already paid out in respect thereof.” Je renvoie à la lecture de Le contrôle de la solvabilité des compagnies d’assurance  en ligne sur le site de l’OCDE, pour une vision globale des approches de ces provisions. La SOA avait publié un rapport en 2009, Comparison of Incurred But Not Reported IBNR Methods que j’encourage à lire.

Nous aborderons mercredi les triangles. Parmi les triangles que nous manipulerons

> source("https://perso.univ-rennes1.fr/arthur.charpentier/bases.R")

qui contient plusieurs fichiers, dont

> 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

ainsi que le triangle évoqué sur http://rworkingparty.wikidot.com/

> OthLiabData = read.csv("http://www.casact.org/research/reserve_data/othliab_pos.csv",header=TRUE, sep=",")
> library(ChainLadder)
> OL = SumData=ddply(OthLiabData,.(AccidentYear,DevelopmentYear,DevelopmentLag),summarise,IncurLoss=sum(IncurLoss_h1-BulkLoss_h1),
+ CumPaidLoss=sum(CumPaidLoss_h1), EarnedPremDIR=sum(EarnedPremDIR_h1))
> LossTri = as.triangle(OL, origin="AccidentYear",
+ dev = "DevelopmentLag", value="IncurLoss")
> Year = as.triangle(OL, origin="AccidentYear",
+ dev = "DevelopmentLag", value="DevelopmentYear")
> TRIANGLE=LossTri
> TRIANGLE[Year>1997]=NA
> TRIANGLE
      dev
origin      1      2      3      4      5      6      7      8      9     10
  1988 128747 195938 241180 283447 297402 308815 314126 317027 319135 319559
  1989 135147 208767 270979 304488 330066 339871 344742 347800 353245     NA
  1990 152400 238665 297495 348826 359413 364865 372436 372163     NA     NA
  1991 151812 266245 357430 400405 423172 442329 460713     NA     NA     NA
  1992 163737 269170 347469 381251 424810 451221     NA     NA     NA     NA
  1993 187756 358573 431410 476674 504667     NA     NA     NA     NA     NA
  1994 210590 351270 486947 581599     NA     NA     NA     NA     NA     NA
  1995 213141 351363 444272     NA     NA     NA     NA     NA     NA     NA
  1996 237162 378987     NA     NA     NA     NA     NA     NA     NA     NA
  1997 220509     NA     NA     NA     NA     NA     NA     NA     NA     NA

Examen intra, éléments de correction

L’énoncé de l’examen intra est en pdf ici et comme annoncé par courriel, la correction de l’intra est dans le pdf en ligne. Comme personne ne semble en désaccord avec les réponses proposées, les notes seront mises en ligne très bientôt. Concertant les questions 18 et 19 quelques compléments d’explications (que je n’avais pas tapé dans le pdf). On avait vu que l’estimateur du maximum de vraisemblance pour une régression de Poisson était asymptotiquement Gaussien,

https://latex.codecogs.com/gif.latex?\widehat{\boldsymbol{\beta}}_{P}\sim\mathcal{N}(\boldsymbol{\beta},V_\infty(\widehat{\boldsymbol{\beta}}_{P}))

(asymptotiquement) avec

https://latex.codecogs.com/gif.latex?V_\infty(\widehat{\boldsymbol{\beta}}_{P})=\left(\sum_{i=1}^n%20\widehat%20Y_i%20\boldsymbol{X}_i\boldsymbol{X}_i%27\right)^{-1}

Quand on a une régression de type binomiale négative, si on note de manière très générale https://latex.codecogs.com/gif.latex?\omega_i=\text{Var}(Y_i|\boldsymbol{X}_i) (on avait vu en cours qu’il existait plusieurs spécifications possibles pour cette variance conditionnelle). Dans ce cas,

https://latex.codecogs.com/gif.latex?\widehat{\boldsymbol{\beta}}_{BN}\sim\mathcal{N}(\boldsymbol{\beta},V_\infty(\widehat{\boldsymbol{\beta}}_{BN}))

avec

https://latex.codecogs.com/gif.latex?V_\infty(\widehat{\boldsymbol{\beta}}_{P})=\left(\sum_{i=1}^n%20\widehat%20Y_i%20\boldsymbol{X}_i\boldsymbol{X}_i%27\right)^{-1}\left[\sum_{i=1}^n%20\omega_i%20\boldsymbol{X}_i\boldsymbol{X}_i\right]\left(\sum_{i=1}^n%20\widehat%20Y_i%20\boldsymbol{X}_i\boldsymbol{X}_i%27\right)^{-1}

Bref, tout dépend fondamentalement de la spécification de la variance conditionnelle. Sous R, c’est la régression binomiale négative de type 1 qui est considérée, i.e.

https://latex.codecogs.com/gif.latex?\omega_i=\text{Var}(Y_i|\boldsymbol{X}_i)=\phi\cdot%20\mathbb{E}(Y_i|\boldsymbol{X}_i)=\phi%20\cdot%20\widehat{Y}_i

On toujours une relation de la forme

https://latex.codecogs.com/gif.latex?\widehat{\boldsymbol{\beta}}_{QP}\sim\mathcal{N}(\boldsymbol{\beta},V_\infty(\widehat{\boldsymbol{\beta}}_{QP}))

avec (en simplifiant un peu)

https://latex.codecogs.com/gif.latex?V_\infty(\widehat{\boldsymbol{\beta}}_{QP})=\phi\cdot\left(\sum_{i=1}^n%20\widehat%20Y_i%20\boldsymbol{X}_i\boldsymbol{X}_i%27\right)^{-1}

aussi, on a

https://latex.codecogs.com/gif.latex?V_\infty(\widehat{\boldsymbol{\beta}}_{QP})=\phi\cdot%20V_\infty(\widehat{\boldsymbol{\beta}}_{P})

Mais comme annoncé en cours, des points étaient données pour ceux qui se contentaient d’affirmer que la variance des estimateurs était plus grande s’il y avait sur-dispersion.

Happy St Patrick’s Day

I love Saint Patrick’s Day for, at least, two reasons. The first one is that, on March 17th, you can play out loud The Pogues, the second one is that it’s the only day in the year when I really enjoy getting a Guiness in a pub. And Guiness is important in statistical science (I did mention a couple of hours ago – on this blog –  that beers were important for social reasons in the academic world, but that was for other reasons…)

> theta=seq(0,pi/2,length=101)
> leaf=sin(2*theta)+.25*sin(6*theta)
> for(k in 0:3)
+ polygon(leaf*cos(theta+k*pi/2),leaf*sin(theta+k*pi/2),col="green")

As mentioned in all my statistics and econometrics courses, the history of statistics (I mean here mathematical statistics) is closely related to Guinness.

A long time ago, there was a Guinness Brewing Company of Dublin, which – as its name suggests – was an Irish brewing company. And the boss, who was to inherit the family business, decided to attract young students, trained in chemistry at Cambridge or Oxford.

In 1899, William Sealy Gosset, who had obtained a double degree in math and chemistry, left Oxford to Dublin. And to be quite honest, being graduate in maths meant when he had studied differential equations and astronomy. Basically, mathematics were useless for Guinness, and he got there with his expertise in chemistry. In fact, William turned out to be also a very good administrator, but this has nothing to do with our story.

William had good memories of his studies in math, and he wondered if he could find a problem to look at. He started studies on workmanship, noting that conditions vary so much (temperature, from hops, malt, manufacturing conditions …) that there were only few consistent data. The “law of errors”  (the central limit theorem) can not apply under these conditions.

In short, Bill (now we know each other a little, we’ll call him Bill) took many measurements, and noticed that the Poisson distribution could be an interesting model to work with. To make the story short, Bill managed to use statistical techniques to control the variance of the production, meaning that he was able to lower losses in the production of beer.

A nice application like this one deserved publication in a scientific journal … Well, of course the Poisson distribution has long been known (it was 1904 and a few months before, Von Bortkiewicz found elegant applications of this law, as discussed in a post  a few weeks ago). But there was a disclosure issue there: Bill’s contract prohibited him from disclosing secrets to the competitors.

Meanwhile, Bill had met Karl Pearson, who was then editor of Biometrika, and encouraged him to publish his results. In 1906, Bill who had helped Guiness to gain a lot of money – doing applied mathematics can be usefull – managed to take a sabbatical to work with Pearson to Galton Laboratory biometrics. Bill and Karl decided to publish the work under a pseudonym “Student.” The legend claims that they had hesitated to use “pupil.”

And for almost 30 years, “Mr Gosset” honorable employee Guinness led a dissolute life by publishing in statistical journals (after work in the brewery) always under the pseudonym “Student”. Of course, it might not be that simple. I mean, Bill had a family life, too. And his wife was the captain of the national Hockey team. So I hardly imagine Bill playing the smart ass and doing mathematical computations, when it was time to wash the dishes or iron his shirt…

In 1908, he wrote a remarkable “the probable error of the mean” remarked, at least, by Ronald Fisher. In fact, Bill found that there was a interesting law, but – as the normal – it was difficult to manipulate to obtain confidence intervals. Without a computer, he had the idea of ​​using monte carlo methods to tabulate quantiles and construct its tables. And he was probably the first one to look carefully at the problem of small samples, unlike Karl Pearson, who always put focus on the asymptotic case.

In fact, looking at his small sample, he saw the denominator magnitudes very close to those specifically manipulated Karl, in particular a square root of chi-square law. Well, of course, remained the normality assumption, but at least we had some results for finite samples !

For the story, William Gosset suggested to use letter z for its statistics, the ratio between the mean and (empirical) standard deviation. But a few years later, statisticians became accustomed to use this letter for Gaussian distribution (i.e. when the variance is known), and it became the standard to use the letter t. Hence finally the present name of “Student-t distribution” and in regression outputs, we have the “t-test”.

A legend (told by Harold Hotelling in his memoirs) claims that the Guinness family discovered this double life on the day of the death of William Gosset in 1937 when mathematicians requested financial assistance to print a volume of the works of their employee. But another legend claims that Mr Guinness himself would have suggested his nickname when he had expressed his intention to publish his research… So I guess we’ll never know. But at least, I’ll think about Bill when I’ll get my first Guiness tonight (but I will probably not be able to tell this story anymore when I’ll reach the fourth…)

Bloguer, et autres activités d’un enseignant-chercheur

J’avais promis un billet bloguer c’est nul dans mon dernier billet sur la recherche et l’activité de blogger. Il va devoir attendre un peu… Dans ce billet, j’avais fait un diagramme pour expliquer la vie et le partage du temps pour un enseignant-chercheur et blogger. Bien entendu, c’est plus compliqué ! Déjà pour un enseignant-chercheur normal, la vie se partage en

  • 50% recherche
  • 50% enseignement
  • 30% tâches administratives diverses

(et oui, c’est dur à représenter dans un diagramme, Tom Roud ou Mix la Malice avaient tenté un décompte du temps passé dans chacune des activités, je renvoie à leurs billets respectifs). Dans ces tâches divers, on retrouve dedans se trouvent les taches officielles, qui figurent sur le CV, comme faire partie du comité bibliothèque ou du comité de programme,voire être responsable d’une programme de maîtrise, mais surtout les tâches non-officielles (en tout cas que parfois on évitera de mettre sur son CV), comme mettre du scotch-tape sur ses factures et discuter avec les services financiers pour comprendre pourquoi le remboursement d’un repas lors d’une conférence est refusé alors que le sandwich à l’aéroport est remboursé ou remplir un sondage en ligne pendant 40 minutes, commandité par un organisme subventionnaire qui veut savoir si le financement de la recherche fonctionne bien (que vous faites car vous voulez expliquer qu’il faut continuer à financer la recherche, mais vous ne trouvez pas la section où vous voulez insister sur le fait que perdre 40 minutes à remplir un sondage ou 2 heures pour vous faire rembourser 6.37 dollars pour un sandwich, c’est un peu du temps pris sur le temps qu’on veut passer à faire de la recherche). Je ne vais pas passer des heures à constituer un bestiaires de ces activités diverses, je renvoie à http://laviedemix.over-blog.com/, le merveilleux blog de mixlamalice, qui fait cela très bien !

Une autre activité (que je compte dans ces tâches diverses) que j’ai pratiqué plusieurs années lorsque j’étais (co)responsable de Master 1 en France, c’est faire jouer les réseaux. Pour faire des recommandations pour des étudiant(e)s qui souhaitent partir dans des programmes plus sélectifs (j’utilise le terme recommandation car il n’y a pas que la lettre que l’on envoie, il y a aussi souvent le coup de fil – ou le courriel – que l’on passe si  on pense vraiment que l’étudiant(e) le mérite), ou pour aider un étudiant à avoir un stage. Et dans ce cas, avoir un réseau, et l’entretenir, ça sert. Pour les stages en entreprises, il y a les copains de promos (j’ai fait une école il y a quelques années), les anciens collègues (j’ai travaillé dans le privé), mais il y a aussi les formations professionnelles que l’on peut faire, les invitations à déjeuner (pour parler boulot, ou à prendre une bière ou deux si on veut en plus le faire sur un mode plus détendu). Dans le monde de la recherche, ce réseau se crée lors de conférences, lors de visites de collègues. Et c’est un peu là que l’on commence à dépasser les 100% de la vie professionnelle dans mon décompte introduisant ce billet: les bons réseaux sont constitués de collègues mais surtout d’ami(e)s. De collègues avec qui ont pourra passer des vacances, qui restent squatter à la maison quand on les invite parler au séminaire (voire faire une semaine de recherche). Ceux qui connaissent l’adresse du blog de mon fils verront dans les photos de vacances des responsables de master, des co-éditeurs de revues, des organisateurs de conférences (voire des moitiés de théorèmes), en maillot de bain ou en train de faire la cuisine… Mais je m’égare. Mon point était que les activité en ligne sont un bon outil pour avoir, ou entretenir, un réseau. Et n’ayant pas de page facebook, le blog est un outil merveilleux pour ça !

Petit exemple pas plus vieux que la semaine dernière. Ma femme (qui est une mathématicienne) recoit un courriel d’un ancien collègue: “j’ai un bon étudiant qui cherche un stage, de recherche, si possible au Canada“. Alors je vais dire un truc avec des pincettes pour éviter de froisser des collègues, mais au Canada, quand on demande un stage à un prof, c’est un peu comme faire la manche de le métro aux heures de pointe à Paris: “salut, j’ai été contactée par un étudiant qui cherche un stage et…“, “oh là là, ça va être dur, j’ai pas l’sous“. En France, ça serait “j’ai pas le temps” mais au  Canada, c’est les sous que l’on mentionne avant tout (cela dit, en France non plus les sous ne sont pas là , mais les règles sont claires et un(e) étudiant(e) qui fait un stage à l’université se doute bien qu’il n’aura rien). Bref, les prises de contact de ma femme n’ont pas été très concluantes… Comme j’évite de parler boulot avec ma femme, il a fallu un peu de temps pour qu’elle évoque cette histoire, en m’expliquant un peu les thèmes que voudraient aborder l’étudiant. Je n’insisterais pas sur le fait que ma première question a été “il a une page internet ?” (ma femme a du me rappeler que les gens normaux n’ont pas de page internet quand ils sont en maîtrise). Mais j’ai glissé “tu sais, il  y a un super prof à Vancouver qui bosse sur ces thèmes là… on se retweete l’un l’autre sur twitter, et il a un super blog… je peux lui envoyer un mot sur twitter, voir si ça l’intéresse…“. Bref, j’envoie un message (privé). Et vendredi, ma femme apprenait que l’étudiant avais pris contact, que le contact semblait bien passé, et donc il devrait être cet été face au Pacifique…

La morale de mon histoire ? L’étudiant avait un CV qui parlait pour lui, et je n’ai pas fait grand chose. Et sans moi, l’étudiant aurait probablement trouvé, malgré tout, un stage passionnant. Mais le blog permet d’entrer dans un réseau de blogueurs (ou de blogueurs-lecteurs de blogs), avec qui on peut discuter, qui permet d’entrer en contact avec des personnes peut-être plus ouvertes que d’autres enseignants-chercheurs. Peut-être est-ce une illusion, peut-être qu’à force de les lire, d’apprendre sur leur vie privée, on a l’illusion de les connaitre… Mais mes expériences ont toujours été bonnes… Les blogs permettent d’entrer en contact avec une communauté de gens passionnants. Les contacts restent souvent virtuels comme lorsque je discute avec David MonniauxBaptiste CoulmontGizmoAlexandre Delaigue, Joel Gombin ou Olivier Bouba-Olga (damned, je me rends compte que je vais en oublier un paquet) mais j’ai toujours du plaisir à aller manger avec Tom Roud, je dois toujours une bière à Reka, et un resto à Camille (entre autres afin de découvrir qui se cache derrière cette identité secrète, c’est mesquin, je sais, mais c’est pas moi qui ait commencé). Derrière les blogueurs se cachent souvent des gens curieux et passionnés (dans le meilleur des cas passionnants), et c’est pour  se rapprocher de cette communauté qu’il est intéressant de bloguer. Bon, j’en conviens, il y a aussi des blogueurs imbus d’eux-même et pénibles… Mais je ne vais pas commencer mon billet bloguer c’est nul aujourd’hui.

Des tablettes… aux blogs ?

 Je viens de finir ‘du scribe au savant‘ d’Yves Gingras, Peter Keating et Camille Limoges. Il m’aura fallu du temps car j’ai l’impression d’avoir de moins en moins de temps pour lire ces temps ci (ou alors des articles ou des rapports qu’il faut évaluer). Pourtant le livre est tout petit quand on y pense. Je veux dire par là qu’on part de la Mésopotamie à l’histoire récente (disons, jusqu’à 1850 environ), en parcourant l’histoire des maths, de la médecine, des sciences naturelles, etc. Forcément, c’est sobre et court, alors quand une section nous intéresse moins, on prend sur soi une page (ou deux maximum), et quand au contraire on se passionne, on reste un peu sur sa faim. J’y ai trouvé une porte ouverte vers plein d’ouvrages même s’il a fallu que je trouve des références par moi-mêmes: à la fin de chaque paragraphe, on attend une référence pour aller plus loin ! C’est un peu ce qui est fait dans le chapitre sur les romains (ce n’est pas le chapitre que j’ai préféré), mais les références sont plus sommaires dans les autres chapitres. J’ai adoré le premier chapitre, qui met en parallèle l’Égypte et de la Mésopotamie, et j’ai vraiment appris plein de choses. J’ai beaucoup aimé l’idée de mettre le scribe (spécialiste de l’écriture et attaché aux inventaires administratifs du royaume) à l’origine de la naissance de la science. Ce qui met un peu l’actuaire à l’origine de la connaissance et du savoir, « actuarius: scribe chargé de la rédaction des procès-verbaux », et j’aime quand on essaye de me raconter une histoire qui me caresse dans le sens du poil. J’ai adoré lire une page (voire deux) sur l’importance de la connaissance des régimes des vents (en fonction des saisons, et des hémisphères) pour les grandes explorations (point dont je n’avais jamais pris conscience, avant de le lire). Ou sur la naissance des universités en Europe (beaucoup moins religieuses que je ne le croyais, mais là encore, à peine deux références…). En plus les schémas sont sobres, et clairs. Ce livre parle du savoir, et de sa transmission (essentiellement en occident, là aussi je suis resté un peu sur ma fin, mai si le livre avait été plus gros, je pense que j’aurais eu peur). J’attends la suite – même si ça ne serait plus réellement un liste d’histoire – qui devrait arriver sur les modes contemporains de transmission du savoir. Avoir un point de vue d’historiens des sciences enrichirait la discussion. Bref, pendant 3 semaines, j’ai appris des choses, tous les jours, et j’avoue qu’on y prend gout !

Examen intra, régression logistique et de Poisson

L’examen intra du cours ACT2040 aura lieu mercredi matin, de 9:00 à 12:00. Aucun document autorisé, sauf les calculatrices (modèle standard, cf plan de cours), et les téléphones seront formellement interdits. Il y aura 34 questions portant sur la première partie du cours (jusqu’à la fin des modèles de comptage, sections 1 à 5 des transparents). 15 questions porteront sur la base décrite dans un précédant billet, sur le nombre de relations extra-conjugales. Il s’agira de décrire les sorties en ligne ici. Je laisse 36 heures pour prendre connaissance de ces sorties. Une version sera donnée lors de l’examen (imprimée 2 pages par feuille, comme dans la version en ligne: si quelqu’un a besoin d’un exemplaire imprimé plus gros, merci de me le faire savoir avant l’examen).

Triangles et provisionnement

La première partie du cours sur le provisionnement (calcul des provisions pour sinistres à payer) aura lieu dans 10 jours. Les transparents sont en ligne ici, et portent sur la construction des triangles de paiements. La méthode chain ladder (et la formalisation proposée par Thomas Mack) ainsi que les extensions seront présentées. La seconde partie portera sur les méthodes basées sur la régression de Poisson.

Génération hackers ?

Je suis régulièrement effaré lorsque j’entends les parents d’ami(e)s de ma fille (la petite dernière) s’extasier que leur enfant arrive à dévérouiller “tout seul !” leur téléphone cellulaire. Le coté auto-célébration de son rejeton, je connais depuis plus d’une dizaine d’années (et je pense en avoir fait aussi ma part), mais c’est plutôt le commentaire “ah vraiment, cette génération sait tellement bien utiliser l’informatique… tu te rends compte, elle a seulement 3 ans !“.

Il serait temps de rétablir une vérité historique: je commence à me faire vieux, mais pourtant, quand j’étais à l’école primaire, on avait des cours d’informatique. Ce n’était pas comme aujourd’hui, où mes enfants ont toujours eu un ordinateur au fond de la salle (dès l’école maternelle, soit à l’age de 2 ans en France), qui servait à faire défiler des photos, ou à aller chercher des informations sur internet. De mon temps, on n’avait pas un ordinateur dans chaque classe, mais il y avait une salle informatique avec des MO5 (qui venaient de sortir). Et par informatique, j’entends par là que l’on apprenait à taper des lignes de code (on n’avait pas internet, oui, j’ai connu ces ages reculés). Ce n’était pas bien méchant: il y avait le logo, mais surtout le Basic (c’était dans le cadre du plan informatique pour tous). De manière surprenante, on était actif devant un écran, en produisant des choses !

Je me souviens avoir fait mes premiers jeux à l’age de mon fils (lui a 10 ans, mais j’étais peut-être un peu plus vieux, maintenant que je découvre les dates exactes en tapant ce billet, je vois que ma mémoire me joue des tours). Souvent, on recopiait des lignes de codes, mais quand même ! Comment se fait-il que l’aspect programmation ait été complètement enlevé des écoles, pour ne parler que d’applications ? et le plus souvent sous une forme ultra rudimentaire, pour ne pas dire mauvaise. Par exemple, je suis régulièrement effaré que mon fils n’ait pas la moindre notion de sécurité informatique (et je ne parle pas de mise en garde avant de créer une page sur Facebook… ce qui serait la moindre des choses) ! Il rentre régulièrement de l’école, m’expliquant que des copains lui ont parlé d’un nouveau site de jeux. Mais il faut s’inscrire, créer un compte, etc. L’autre jour, en discutant un peu avec lui des mots de passe, on a compris qu’il pourrait être amusant de créer un site, où les gens entreraient leur adresse électronique, et un mot de passe. Et que le mot de passe permettrait surement de se connecter sur leur adresse électronique (oui, il s’est rendu compte qu’il utilisait toujours le même mot de passe). Il pourrait alors faire des bonnes blagues en envoyant des messages à leur place ! (il a 10 ans, ses blagues sont encore relativement innocentes). On a eu une longue discussion au début de la semaine de relâche, et je me rends compte qu’il ne manque pas grand chose pour en faire un vrai hacker (j’entends par là le sens que l’on trouve sur wikipedia, “un hacker est quelqu’un qui aime comprendre le fonctionnement d’un mécanisme, afin de pouvoir le bidouiller“; il est dommage que dans l’imaginaire collectif – j’entends par là ce que l’on peut lire dans les journaux – le hacker soit aussi mal perçu, alors que c’est juste quelqu’un de curieux… la curiosité est devenu un bien mauvais défaut).

Après avoir longuement hésité (et parce que je ne trouvais pas de camps de jour qui proposait d’apprendre à faire de l’informatique), je me suis lancé: mon fils va faire du R. Lui est content car il réclame souvent à pouvoir “faire des trucs” sur l’ordinateur, et moi car j’ai l’illusion qu’il va apprendre des choses qui pourraient lui servir un jour (au moins à comprendre comment fonctionne un ordinateur). Et autant que faire se peut, j’essaye de séparer les activités familiales de ce qui pourrait s’approcher du travail. En fait, je pensais acheter python for kids (et en profiter pour découvrir un langage que j’aurais du apprendre voilà 10 ans). Mais le livre est en anglais, et mon fils n’est pas très à l’aise. Bref, je me suis lancé dans R for kids (par moi même)L’objectif était d’apprendre à faire des dessins  (un peu dans l’idée du logo je pense). De comprendre qu’un dessin était une succession de formes de base. J’ai commencé (pendant le cours d’escrime de mon fils) à taper quelques fonctions simples (carré, trait, triangle, disque, point, etc), et à coder les principales couleurs (pour qu’il les tape en français). Tout est caché dans la fonction

source("http://freakonometrics.free.fr/RforKIDZ.R")

Ensuite, on s’est lancé. Le point de départ est de faire un dessin ! Oui, c’est plus simple. Ensuite, on va définir des points, en donnant leur coordonnées à partir de la grille de fond

fond()

pour créer la grille de fond, et pour les points

A=c(0,0)
B=c(4,12)
C=c(0,12)
D=c(2,15)

On peut d’ailleurs visualiser ces points (pour vérifier qu’ils sont bien placés)

point(A)
point(B)
point(C)
point(D)

Puis on fait les figures.

carre(A,B,"gris")
polygone(C,B,D,couleur="rouge")

(dans la première version, je n’avais pas pensé faire une fonction spécifique pour les triangles) Ensuite, on rajoute un drapeau,

E=c(2,18)

pour le sommet du mat, puis pour le reste

trait(D,E,"noir",ep=2)
F=c(2,16)
point(F)
G=c(6,17)
point(G)
polygone(E,F,G,couleur="jaune")

Ce n’est pas du code, ça…. Ben, un peu quand même…. surtout quand on a vu qu’on pouvait translater une figure,

h=15
A=c(0+h,0)
B=c(4+h,12)
C=c(0+h,12)
D=c(2+h,15)
E=c(2+h,18)
carre(A,B,"gris")
polygone(C,B,D,couleur="rouge")
trait(D,E,"noir",ep=2)
F=c(2+h,16)
G=c(6+h,17)
polygone(E,F,G,couleur="jaune")

Et hop, on a deux tours.

Je pense que c’est là le cœur de la programmation: comprendre qu’il y a une forme de base, et qu’après, on ne fait que répéter. Ensuite, au centre, on a fait le mur, et on a mis les créneaux, là encore, en comprenant que c’était la même figure, translatée plusieurs fois…

A=c(4,0)
B=c(15,8)
carre(A,B,"gris")

A=c(5,8)
B=A+1
carre(A,B,"gris")

et puis on répète

A=c(7,8)
B=A+1
carre(A,B,"gris")

A=c(9,8)
B=A+1
carre(A,B,"gris")

etc…

(on a vu au passage que si on ne compte pas utiliser un point, on peut donner son nom à un autre) Enfin, pour faire le mur et la porte, on a vu qu’on avait le choix: Le plus simple (après de longues négociations discussions) a été de faire un rectangle, puis de faire un trou carré, et un cercle (en blancs).

A=c(13,8)
B=A+1
carre(A,B,"gris")

A=c(8,0)
B=c(12,4)
carre(A,B,"blanc")

C=c(10,4)
disque(C,2,"blanc")

Pas mal ? Bon, maintenant, ça me gêne un peu. Parce que je ne suis pas un bon codeur, et je vais apprendre de mauvais réflexes à mon fils (ou mes enfants, car ma fille a fini par participer, mais on verra plus tard ce qu’elle a fait).

Le plus marrant, c’est qu’on a vu comment faire un film: on a construit une voiture (assez sommaire, j’en conviens, les ingénieurs pesteront surement en voyant notre boite à chaussures avec deux roues).

dessin=function(x){
fond()
A=c(2+x,2)
B=c(6+x,4)
C=c(3+x,2)
D=c(5+x,2)
carre(A,B,"vert")
disque(C,.75,"noir")
disque(D,.75,"noir")
}

(c’est moi qui est codé la fonction, on verra ça plus en détails une prochaine fois). Et ensuite, on l’a faite se translater, de la gauche, vers la droite: on commence par taper

dessin(0)

puis

dessin(1)

et

dessin(2)

etc… en allant vite, on crée du mouvement…

Amusant, non ?

Ma fille, elle a opté pour un dessin plus traditionnel… le fameux “maison avec arbre et arc en ciel”,

Voilà ce qu’on a pu faire en quelques heures… Je pense qu’on pourrait faire mieux, et je suis ouvert à toutes suggestions: sur la façon d’apprendre à coder, sur l’interface (on utilise R-studio: on code dans la fenêtre de gauche, on utilise l’icône “run” et le dessin s’affiche à droite), sur d’éventuelles applications amusantes, ou des expériences menées par des instituteurs qui veulent apprendre les bases de l’informatique à leurs élèves. Je trouve énormément de ressources en anglais, comme le livre python for kids dont je parlais au début (ou les sites dédiés aux jeux que l’on code soi-même, en python, comme inventwithpython.com/, qui me font penser à ce que je faisais quand j’étais petit), mais je pourrais citer le scratch. Car si la France semble avoir été pionnière en 1980, je ne vois plus grand chose ces temps-ci, en français… Mais je ne sais peut-être pas bien chercher dans la communauté francophone.

Multiple (smoothed) regression and portfolio exposure

Wednesday, in class, we’ve seen how to visualize a multiple regression model (with two continuous explanatory variables). Here, the goal is to predict the average cost of an insurance claim, using some covariates, e.g. the age of the driver, and the age of the car (recall that losses here are liability losses). The prediction obtained from a (standard) generalized linear model, with a log-link

> reg1=glm(cout~ageconducteur+agevehicule,data=base,family=Gamma(link="log"))

The code to visualize the predicted average cost is the following: first, we have to compute predictions for specific values,

> pred=function(x,y){
+ predict(reg,newdata=data.frame(ageconducteur=x,
+ agevehicule=y),type="response")

Then, we use this function to compute values on a grid,

> X=seq(20,80,by=5)
> Y=0:20
> Z=outer(X,Y,p)
> image(X,Y,Z,col=rev(heat.colors(101)))
> contour(X,Y,Z,add=TRUE,
+ levels=c(1400,1800,2000,2200,2400,2600,2800,3000,3200,4000,5000))

If we use factors, and not continuous variates (cut versions of those two variates),

> reg2=glm(cout~cut(ageconducteur,breaks=c(0,22,35,55,80,100))*
+               cut(agevehicule,breaks=c(-1,1,3,5,10,100)),
+ data=base,family=Gamma(link="log"))

(note that we consider the Cartesian product, so values are computed for each product of factors, age of the driver and age of the car) we obtain

Obviously, we’re missing something here: the most expensive class with one model is the cheapeast for the other one! Of course, it might come from our classes (that were chosen a bit randomly), but it might be interesting to use nonlinear functions of the ages. So, let us use splines to smooth those two variables,

> reg3=glm(cout~bs(ageconducteur)+bs(agevehicule),data=base,
+ family=Gamma(link="log"))

With additive smoothed functions, we obtained a symmetric graph (due to the additive property)

while with a bivariate spline

> library(mgcv)
+ reg4=gam(cout~s(ageconducteur,agevehicule),data=base,
+ family=Gamma(link="log"))

(for some odd reasons, I could not use – easily – a bivariate spline in the Generalized Linear Model, but it did work considering a Generalized Additive Model – which is, by no means additive now). We can identify here some regions where the average cost can be extremely expensive… But, as mentioned wednesday, one should keep in mind that some parts of the square above are not reached. More precisely, the distribution of the portfolio, as a function of those two covariates is the following

Thus, the proportion of young drivers driving a brand new car, and the proportion of old drivers driving a very old car is rather small… If the goal is to find niches, one should look at the prediction more carefully, but if the goal is to make that everyone gets an insurance cover, maybe we should allow that some drivers are under-priced (especially when are rare in the portfolio). And one should keep in mind that average costs are extremely sensitive to large losses, as discussed previously http://freakonometrics.hypotheses.org/3490 (and in class)

In the univariate case, I have migrated an old post, we I tried to reproduce (in R and in French) some standard graphs in the insurance industry: it is always interesting to visualize not only the prediction obtained from our models, but also the size of each class in the portfolio,

The post is online here http://freakonometrics.hypotheses.org/1224

Comparing quantiles for two samples

Recently, for a research paper, I got some samples, and I wanted to compare them. Not to compare their means (by construction, all of them were centered) but there dispersion. And not their variance, but more their quantiles. Consider the following boxplot type function, where everything here is quantile related (which is not the case for standard boxplot, see http://freakonometrics.hypotheses.org/4138, in French)

> boxplotqbased=function(x){
+ q=quantile(x[is.na(x)==FALSE],c(.05,.25,.5,.75,.95))
+ plot(1,1,col="white",axes=FALSE,xlab="",ylab="",
+ xlim=range(X),ylim=c(1-.6,1+.6))
+ polygon(c(q[2],q[2],q[4],q[4]),1+c(-.4,.4,.4,-.4))
+ segments(q[1],1-.4,q[1],1+.4)
+ segments(q[5],1,q[4],1)
+ segments(q[5],1-.4,q[5],1+.4)
+ segments(q[1],1,q[2],1)
+ segments(q[3],1-.4,q[3],1+.4,lwd=2)
+ xt=x[(x<q[1])|(x>q[5])]
+ points(xt,rep(1,length(xt)))
+ axis(1)
+ }

(one can easily adapt the code for lists, e.g.). Consider for instance temperature, when the (linear) trend is removed (see http://freakonometrics.hypotheses.org/1016 for a discussion on that series, in Paris),

from January 1st till December 31st. Let us remove now the seasonal cycle, i.e. we do have here the difference with the average seasonal temperature (with here upper and lower quantiles),

Seasonal boxplots are here (with Autumn on top, then Summer, Spring and Winter, below),

If we zoom in, we do have (where upper and lower segments are 95% and 5% quantiles, while classically, boxes are related to the 75% and the 25% quantiles)

Is there a (standard) test to compare quantiles – some of them perhaps ? Can we compare easily quantiles when have two (or more) samples ?

Note that this example on temperature could be related to other old posts (see e.g. http://freakonometrics.hypotheses.org/2190), but the research paper was on a very different topic.

Consider two (i.i.d.) samples https://latex.codecogs.com/gif.latex?\{x_1,\cdots,x_m\} and https://latex.codecogs.com/gif.latex?\{y_1,\cdots,y_n\}, considered as realizations of random variables https://latex.codecogs.com/gif.latex?X and https://latex.codecogs.com/gif.latex?Y. In all statistical courses, tests on the average are always considered, i.e.

https://latex.codecogs.com/gif.latex?H_0:\mathbb{E}(X)=\mathbb{E}(Y)

against

https://latex.codecogs.com/gif.latex?H_1:\mathbb{E}(X)\neq\mathbb{E}(Y)

Usually, the idea in courses is to start with a one sample test, and to test something like

https://latex.codecogs.com/gif.latex?H_0:\mathbb{E}(X)=\mu_\star

against

https://latex.codecogs.com/gif.latex?H_1:\mathbb{E}(X)\neq\mu_\star

The idea is to assume that samples are from Gaussian variables,

https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20T%20=%20\frac{\overline{x}%20-%20\mu_\star}{\widehat{\sigma}/\sqrt{n}}
Under https://latex.codecogs.com/gif.latex?H_0https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20T has a Student t distribution. All that can be found in any Statistics 101 course. We can derive https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20pvalue, computing probabilities that https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20T exceeds the observed values (for two sided tests, the probability that the absolute value of https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20T exceed the absolute value of the observed statistics). This test is closely related to the construction of confidence intervals for https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20\mu. If https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20\mu_\star belongs to the confidence interval, then it might be a suitable value. The graphical representation of this test is related to the following graph

Here the observed value was 1,96, i.e. the https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20pvalue (the area in red above) is exactly 5%.

To compare means, the standard test is based on

https://latex.codecogs.com/gif.latex?T%20=%20{\overline{x}%20-%20\overline{y}%20\over%20%20\displaystyle\sqrt{{s_x^2%20\over%20m}%20+%20{s_y^2%20\over%20n}}%20}

which has – under https://latex.codecogs.com/gif.latex?H_0 – a Student-t distribution, with https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20\nu degrees of freedom, where

https://latex.codecogs.com/gif.latex?\nu%20=%20\frac{(s_x^2/m%20+%20s_y^2/n)^2}{(s_x^2/m)^2/(m-1)%20+%20(s_y^2/n)^2/(n-1)}.

Here, the graphical representation is the following,

But tests on quantiles are rarely considered in statistical courses. In a general setting,define quantiles as

https://latex.codecogs.com/gif.latex?Q_X(p)=\inf\left\{%20x\in%20\mathbb%20R%20:%20p%20\le%20\mathbb%20P(X\leq%20x)%20\right\}

one might be interested to test
https://latex.codecogs.com/gif.latex?H_0:Q_X(p)=Q_Y(p)
against
https://latex.codecogs.com/gif.latex?H_1:Q_X(p)\neq%20Q_Y(p)
for some https://latex.codecogs.com/gif.latex?p\in(0,1). Note that we might be interested also to test if

https://latex.codecogs.com/gif.latex?H_0:Q_X(p_k)=%20Q_Y(p_k)
for all https://latex.codecogs.com/gif.latex%20?%20%20k, for some vector of probabilities https://latex.codecogs.com/gif.latex?\boldsymbol{p}=(p_1,\cdots,p_d)\in(0,1)^d.
One can imagine that this multiple test will be more complex. But more interesting, e.g. a test on boxplots (are the four quantiles equal ?).  Let us start with something a bit more simple: a test on quantiles for one sameple, and the derivation of a confidence interval for quantiles.

  • Quantiles for one sample

The important idea here is that it should be extremely simple to get https://latex.codecogs.com/gif.latex?pvalues. Consider the following sample, and let us run a test to assess if the median can be zero.

> set.seed(1)
> X=rnorm(20)
> sort(X)
[1] -2.21469989 -0.83562861 -0.82046838 -0.62645381 -0.62124058 -0.30538839
[7] -0.04493361 -0.01619026  0.18364332  0.32950777  0.38984324  0.48742905
[13]  0.57578135  0.59390132  0.73832471  0.82122120  0.94383621  1.12493092
[19]  1.51178117  1.59528080
> sum(X<=0)
[1] 8

Here, 8 observations (out of 20, i.e. 40%) were below zero. But we do know the distribution of https://latex.codecogs.com/gif.latex%20?%20%20N the number of observation below the target

https://latex.codecogs.com/gif.latex?N=\sum_{i=1}^n%20\boldsymbol{1}(X_i\leq%20x_\star)

It is a binomial distribution. Under https://latex.codecogs.com/gif.latex?H_0, it is a binomial distribution https://latex.codecogs.com/gif.latex?\mathcal{B}(n,p_\star) where https://latex.codecogs.com/gif.latex?p_\star is the probability target (here 50% since the test is on the median). Thus, one can easily compute the https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20p-value,

> plot(n,dbinom(n,size=20,prob=0.50),type="s",xlab="",ylab="",col="white")
> abline(v=sum(X<=0),col="red")
> for(i in 1:sum(X<=0)){
+ polygon(c(n[i],n[i],n[i+1],n[i+1]),
+ c(0,rep(dbinom(n[i],size=20,prob=0.50),2),0),col="red",border=NA)
+ polygon(21-c(n[i],n[i],n[i+1],n[i+1]),
+ c(0,rep(dbinom(n[i],size=20,prob=0.50),2),0),col="red",border=NA)
+ }
> lines(n,dbinom(n,size=20,prob=0.50),type="s")

which yields

Here, the https://latex.codecogs.com/gif.latex%20?%20%20%20%20%20p-value is

> 2*pbinom(sum(X<=0),20,.5)
[1] 0.5034447

Here the probability is easy to compute. But one can observe that there is some kind of disymmetry here. Actually, if the observed value was not 8, but 12, some minor changes should be done (to keep some symmetry),

> plot(n,dbinom(n,size=20,prob=0.50),type="s",xlab="",ylab="",col="grey")
> abline(v=20-sum(X<=0),col="red")
> for(i in 1:sum(X<=0)){
+ polygon(c(n[i],n[i],n[i+1],n[i+1])-1,
+ c(0,rep(dbinom(n[i],size=20,prob=0.50),2),0),col="red",border=NA)
+ polygon(21-c(n[i],n[i],n[i+1],n[i+1])-1,
+ c(0,rep(dbinom(n[i],size=20,prob=0.50),2),0),col="red",border=NA)
+ }
> lines(n-1,dbinom(n,size=20,prob=0.50),type="s")

Based on those observations, one can easily write a code to test if the https://latex.codecogs.com/gif.latex?p_\star-quantile of a sample is https://latex.codecogs.com/gif.latex?x_\star. Or not. For a two sided test, consider

> quantile.test=function(x,xstar=0,pstar=.5){
+ n=length(x)
+ T1=sum(x<=xstar)
+ T2=sum(x< xstar)
+ p.value=2*min(1-pbinom(T2-1,n,pstar),pbinom(T1,n,pstar))
+ return(p.value)}

Here, we have

> quantile.test(X)
[1] 0.5034447

Now, based on that idea, due to the duality between confidence intervals and tests, one can easily write a function that computes confidence interval for quantiles,

> quantile.interval=function(x,pstar=.5,conf.level=.95){
+ n=length(x)
+ alpha=1-conf.level
+ r=qbinom(alpha/2,n,pstar)
+ alpha1=pbinom(r-1,n,pstar)
+ s=qbinom(1-alpha/2,n,pstar)+1
+ alpha2=1-pbinom(s-1,n,pstar)
+ c.lower=sort(x)[r]
+ c.upper=sort(x)[s]
+ conf.level=1-alpha1-alpha2
+ return(list(interval=c(c.lower,c.upper),confidence=conf.level))}
> quantile.interval(X,.50,.95)
$interval
[1] -0.3053884  0.7383247

$confidence
[1] 0.9586105

Because of the use of non-asymptotic distributions, we can not get exactly a 95% confidence interval. But it is not that bad, here.

  • Comparing quantiles for two samples

Now, to compare quantiles for two samples… it is more complicated. Exact tests are discussed in Kosorok (1999) (see http://bios.unc.edu/~kosorok/…) or in Li, Tiwari and Wells (1996) (see http://jstor.org/…). For the computational aspects, as mentioned in a post published almost one year ago on http://nicebread.de/… there is a function to compare quantiles for two samples.

> install.packages("WRS")
> library("WRS")

Some multiple tests on quantiles can be performed here. For instance, on the temperature, if we compare quantiles for Winter and Summer (on only 1,000 observations since it can be long to run that function), i.e. 5%, 25%, 75% and 95%,

> qcomhd(Z1[1:1000],Z2[1:1000],q=c(.05,.25,.75,.95))
q   n1   n2      est.1      est.2 est.1_minus_est.2     ci.low     ci.up     p_crit p.value signif
1 0.05 1000 1000 -6.9414084 -6.3312131       -0.61019530 -1.6061097 0.3599339 0.01250000   0.220     NO
2 0.25 1000 1000 -3.3893867 -3.1629541       -0.22643261 -0.6123292 0.2085305 0.01666667   0.322     NO
3 0.75 1000 1000  0.5832394  0.7324498       -0.14921041 -0.4606231 0.1689775 0.02500000   0.338     NO
4 0.95 1000 1000  3.7026388  3.6669997        0.03563914 -0.5078507 0.6067754 0.05000000   0.881     NO

or if we compare quantiles for Winter and Summer

> qcomhd(Z1[1:1000],Z3[1:1000],q=c(.05,.25,.75,.95))
q   n1  n2      est.1     est.2 est.1_minus_est.2     ci.low       ci.up     p_crit p.value signif
1 0.05 1000 984 -6.9414084 -6.438318        -0.5030906 -1.3748624  0.39391035 0.02500000   0.278     NO
2 0.25 1000 984 -3.3893867 -3.073818        -0.3155683 -0.7359727  0.06766466 0.01666667   0.103     NO
3 0.75 1000 984  0.5832394  1.010454        -0.4272150 -0.7222362 -0.11997409 0.01250000   0.012    YES
4 0.95 1000 984  3.7026388  3.873347        -0.1707078 -0.7726564  0.37160846 0.05000000   0.539     NO

(the following graphs are then plotted)

Those tests are based on the procedure proposed in Wilcox, Erceg-Hurn,  Clark and Carlson (2013), online on http://tandfonline.com/…. They rely on the use of bootstrap samples. The idea is quite simple actually (even if, in the paper, they use Harrell–Davis estimator to estimate quantiles, i.e. a weighted sum of ordered statistics – as described in http://freakonometrics.hypotheses.org/1755 – but the idea can be understood with any estimator): we generate several bootstrap samples, and compute the median for all of them (since our interest was initially on the median)

>  Q=rep(NA,10000)
>  for(b in 1:10000){
+  Q[b]=quantile(sample(X,size=20,replace=TRUE),.50)
+  }

Then, to derive a confidence interval (with, say, 95% confidence), we compute quantiles of those median estimates,

> quantile(Q,c(.025,.975))
     2.5%     97.5% 
-0.175161  0.666113

We can actually visualize the distribution of that bootstrap median,

> hist(Q)

Now, if we want to compare medians from two independent samples, the strategy is rather similar: we bootstrap the two samples – independently – then compute the median, and keep in mind the difference. Then, we will look if the difference is significantly different from 0. E.g.

> set.seed(2)
> Y=rnorm(50,.6)
> QX=QY=D=rep(NA,10000)
> for(b in 1:10000){
+ QX[b]=quantile(sample(X,size=length(X),replace=TRUE),.50)
+ QY[b]=quantile(sample(Y,size=length(Y),replace=TRUE),.50)
+ D[b]=QY[b]-QX[b]
+ }

The 95% confidence interval obtained from the bootstrap difference is

> quantile(D,c(.025,.975))
      2.5%      97.5% 
-0.2248471  0.9204888

which is rather close to was can be obtained with the R function

> qcomhd(X,Y,q=.5)
    q n1 n2    est.1     est.2 est.1_minus_est.2    ci.low     ci.up p_crit p.value signif
1 0.5 20 50 0.318022 0.5958735        -0.2778515 -0.923871 0.1843839   0.05    0.27     NO

(where the difference is here the oppositive of mine). And when testing for 2 (or more) quantiles, Bonferroni method can be used to take into account that those tests cannot be considered as independent.