Category Archives: Tests

Detecting distributions with infinite mean

In a post I published a few month ago (in French, here, based on some old paper,there), I mentioned a statistical procedure to test if the underlying distribution https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO02.pngof an i.i.d. sample https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO01.png had a finite mean (based on extreme value results). Since I just used it on a small dataset (yes, with real data), I decided to post the R code, since it is rather simple to use. But instead of working on that dataset, let us see what happens on simulated samples. Consider https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO03.png=200 observations generated from a Pareto distribution

https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO04.png

(here https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=2, as a start)

> a=2
> X=runif(200)^(-1/a)

Here, we will use the package developped by Mathieu Ribatet,

> library(RFA)
  • Using Generalized Pareto Distribution (and LR test)

A first idea is to fit a GPD distribution on observations that exceed a threshold https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO11.png>1.
Since we would like to test https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO06.png (against the assumption that the expected value is finite, i.e. https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO07.png), it is natural to consider the likelihood ratio

https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO08.png

Under the null hypothesis, the distribution of https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO09.png should be chi square distribution with one degree of freedom. As mentioned here, the significance level is attained with a higher accuracy by employing Bartlett correction (there). But let  us make it as simple as possible for the blog, and use the chi-square distribution to derive the p-value.
Since it is rather difficult to select an appropriate threshold, it can be natural (as in Hill’s estimator) to consider https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO10.png, and thus, to fit a GPD on the https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO13.pnglargest values. And then to plot all that on a graph (like the Hill plot)

> Xs=rev(sort(X))
> s=0;G=rep(NA,length(Xs)-14);Gsd=G;LR=G;pLR=G
> for(i in length(X):15){
+ s=s+1
+ FG=fitgpd(X,Xs[i],method="mle")
+ FGc=fitgpd(X,Xs[i],method="mle",shape=1)
+ G[s]=FG$estimate[2]
+ Gsd[s]=FG$std.err[2]
+ FGc$fixed
+ LR[s]=FGc$deviance-FG$deviance
+ pLR[s]=1-pchisq(LR[s],df=1)
+ }

Here we keep the estimated value of the tail index, and the associated standard deviation of the estimator, to draw some confidence interval (assuming that the maximum likelihood estimate is Gaussian, which is correct only when n is extremely large). We also calculate the deviance of the model, the deviance of the constrained model (https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO06.png), and the difference, which is the log likelihood ratio. Then we calculate the p-value (since under https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO12.png the likelihood ratio statistics has a chi-square distribution).
If https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=2, we have the following graph, with on top the p-value (which is almost null here), the estimation of the tail index he https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO13.png largest values (and a confidence interval for the estimator),

If https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1.5 (finite mean, but infinite variance), we have

i.e. for those two models, we clearly reject the assumption of infinite mean (even if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png gets too close from 1, we should consider thresholds large enough). On the other hand, if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1 (i.e. infinite mean) we clearly accept the assumption of infinite mean (whatever the threshold),

  • Using Hill’s estimator

An alternative could be to use Hill’s estimator (with Alexander McNeil’s package). See here for more details on that estimator. The test is simply based on the confidence interval derived from the (asymptotic) normal distribution of Hill’s estimator,

> library(evir)
> Xs=rev(sort(X))
> HILL=hill(X)
> G=rev(HILL$y)
> Gsd=rev(G/sqrt(HILL$x))
> pLR=1-pnorm(rep(1,length(G)),mean=G,sd=Gsd)

Again, if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=2, we clearly rejct the assumption of infinite mean,

and similarly, if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1.5 (finite mean, but infinite variance)

Here the test is more robust than the one based on the GPD. And if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=1 (i.e. infinite mean), again we accept https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO12.png,

Note that if https://perso.univ-rennes1.fr/arthur.charpentier/latex/paretoRO05.png=0.7, it is still possible to run the procedure, and hopefully, it suggests that the underlying distribution has infinite mean,

(here without any doubt). Now you need to wait a few days to see some practical applications of the idea (there was on in the paper mentioned above actually, on business interruption insurance losses).

Peut on faire l’économie du formalisme quand on parle d’extrêmes ?

Tous les blogs économiques saluent la parution en poche du joli petit livre de Daniel Zajdenweber, Economie des Extrêmes. En particulier, beaucoup de monde salue ce livre qui explique simplement des choses complexes…. Par exemple Alexandre dès 2001 “Passé le premier chapitre, un peu ardu, et qui nécessite du lecteur des connaissances de base en statistique et probabilités (notion de lois de probabilité, d’espérance, de variance…) qui décrit en termes littéraires les caractéristiques de ces lois, l’auteur applique ces résultats à un grand nombre de phénomènes concrets, et en tire les conséquences“. Mais peut-on parler d’économie des extrêmes sans être technique ?

Histoire que mon message ne soit pas déformé, je trouve passionnant ce petit livre introductif à la problématique des risques extrêmes (qui est un de mes dadas depuis quelques années) mais j’espère qu’il servira d’encouragement à une lecture d’ouvrages plus détaillés sur le sujet. Car la vulgarisation a des limites que l’on atteint vite quand on parle de sujets aussi complexes.

L’exemple que j’ai le plus étudié est celui des sinistres de perte d’exploitation (longuement évoqué par Daniel Zajdenweber dans son livre). Il y a quelques années j’avais utilisé cette partie du livre comme base pour faire un sujet d’examen pour le cours de “réassurance et grands risques” que je donnais alors à l’ENSAE1. Et malheureusement, mes compétences littéraires sont très limitées, donc je vais faire des maths. Dans le livre, la figure suivante est présentée,

qui correspond effectivement à la fonction tracée dès 1925 par Karl Gustav Hagstroem (j’avais souligné (ici) ses travaux précurseurs où l’intérêt de la loi de Pareto pour modéliser les très grands riques apparaissait pour la première fois). C’est en effet assez naturel: si on a une loi de Pareto, i.e.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z01.png

alors on pourrait écrire, en passant au logarithme

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z02.png

Si on représente la version empirique, c’est à dire le nuage de points

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z03.png

alors pour une loi de Pareto, les points devraient être alignés suivant une droite, et la pente doit correspondre au paramètre de la fonction puissance. C’est visiblement l’idée exploitée ici.

Autrement dit, les pointillés ne sont un intervalle de confiance, mais juste un outils graphique pour se demander si la pente vaut 1, ou pas. Daniel Zajdenweber affirme que la pente doit ici être -1.

Le fait que la valeur soit unitaire ou pas a en effet un impact très important en terme d’assurabilité du risque de perte d’exploitation. Rappelons que pour une variable positive (et c’est le cas ici). Et si on a une telle loi de Pareto (de puissance unitaire), alors la prime pure d’un traité de réassurance, couvrant entre m et M s’écrit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z06.png

soit

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z07.png

ce qui correspond aux calculs de Daniel Zajdenweber… Mais encore une fois “l’absence d’espérance mathématique de la distribution des sinistres” est une conclusion très forte sur laquelle on peut essayer de revenir.

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z04.png

aussi ici

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z05.png

autrement dit, l’espérance est finie si la pente est strictement plus grande (en valeur absolue) que 1. Si la pente est inférieure (ou égale) à 1, le risque n’est pas assurable ! Ce qui est une conclusion très très forte pour les assureurs.
J’ai donc demandé à la FFSA la base de données utilisée ici, et pour éviter des problèmes d’inflation des coûts de sinistres entre 1992 et 2000. Si je prends tous les sinsitres, on obtient l’ajustement de Pareto suivant

soit une pente (en valeur absolue de 1.47). Mais encore une fois, l’ajustement de Pareto se fait sur les grands sinistres. Hill a proposé un estimateur très populaire pour estimer ce coefficient, où on ne prend en compte que les k observations les plus grandes, et on regarde l’estimation de la pente du graphique de Pareto pour ces quelques valeurs. On représente alors l’estimation en fonction du nombre de grands sinistres, ou du seuil définissant les graphs sinistres. Numériquement, en posant

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z08.png

on peut écrire comme estimateur de la pente

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z09.png

soit, en simplifiant le numérateur,

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z10.png

tel que l’a construit Hill en 1975. Graphiquement, on a ici

Bref, la question est ce savoir si on atteint la valeur 1 pour les grands sinistres. Graphiquement, on a malgré tout envie de rejeter cette hypothèse.
Une solution peut être de faire un test statistique, basé sur de ratio de vraisemblance, comme le suggèrent Reiss & Thomas (2001) or Coles (2001). En fait, on peut même utiliser d’autres estimateurs que celui de Hill, comme celui obtenu en faisant un ajustement de loi GPD (Pareto généralisée
) sur la loi des Excès, ou une loi GEV sur des maximas par blocs (Generalized Extreme Value). On introduit alors la statistique de test suivante

https://perso.univ-rennes1.fr/arthur.charpentier/latex/z12.png

et on regarde les p-value (ainsi que la correction de Bartlett à droite),

On peut aussi, plus simplement, estimer plusieurs coefficients de pentes pour des seuils différents, et regarder la borne supérieure de l’intervalle de confiance,

Bref, même si avec un des ajustements de loi GPD on hésite à retenir une pente unitaire, la plupart des tests rejettent cette hypothèse, et donc le risque de perte d’exploitation semble être assurable, d’espérance mathématique finie. Bref, les dessins c’est très bien pour faire passer une idée, mais ne retenir que ça pour en tirer des conclusions aussi fortes me laisser sceptique….

Choix multiples et théorie des tests

En poursuivant la première question du QCM de “probabilité & statistiques” (en ligne ici pour l’énoncé, et pour une ébauche de correction), on a l’occasion de voir une application intéressante (ou tout du moins selon mon point de vue) de la théorie des tests.

Je souhaite tester une hypothèse

H0:  la personne a répondu au hasard aux questions

(je spécifierais plus tard l’hypothèse alternative, de toutes les façons, on n’en a pas besoin pour construire la statistique de test.)
En répondant au hasard (i.e. sous l’hypothèse H0), la personne avait 1 chance sur 3 de gagner 1 point, et 2 chances sur 3 de perdre 0,5 point1. Bref, à une transformation affine prêt, pour chaque question, on a une loi de Bernoulli. Plus précisément le nombre de points gagnés (ou perdus) à la question i est

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-1.png

Aussi, la loi du nombre total de points est

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-2.png

or si l’on suppose que l’on réponde au hasard à toutes les questions, on peut suppose les variables indépendantes, et donc

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-3.png

et sans être trop prétentieux, on connaît de cette loi (en particulier sa fonction de répartition et sa fonction quantile, au moins numériquement en tous les cas). Bref, si F est cette fonction de répartition,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-6.png

ensuite on sait faire: c’est juste trouver la loi d’une transformation croissante de cette variable (linéaire en fait). Aussi

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-7.png

c’est à dire

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-9.png

Ben ça peut se représenter facilement avec R
> X=seq(-10,20,by=.1)
> G = function(x){pbinom((10+x)*2/3,size=20,prob=1/3)}
> Y=G(X)
> plot(X,Y,type=”l”)

Notons que je peux aussi avoir la fonction quantile de cette loi. On cherche q tel que

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-10.png

C’est à dire tel que

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-11.png

i.e.

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-12.png

En  R, on obtient la fonction suivante
> G.inv = function(u){1.5*qbinom(u,size=20,prob=1/3)-10}
Voilà, j’ai donc la loi de la note (qui me servira de statistique de test) sous l’hypothèse H0. Bon, ça suffit presque à construire complètement le test. Il suffit de spécifier l’hypothèse alternative afin de construire les régions critiques. A priori je vois deux formes pour l’hypothèse alternative,

H1:  la personne n’a pas répondu au hasard aux questions

ce qui correspondra à un test bilatéral, incluant ceux qui ont des bonnes notes (et qui sont plus fort que le hasard), et ceux qui ont des très mauvaises notes (et qui sont donc moins fort que le hasard en quelque sorte). L’autre hypothèse alternative, qui a plus de sens à mes yeux serait

H1:  la personne est plus forte que le hasard

ce qui sous-entend qu’elle a un peu travaillé (ou qu’elle a raisonné un minimum).
Dans le premier cas, on cherche les deux seuils de la région de rejet
> (note.inf = G.inv(0.05))
-5.5
> (note.sup=G.inv(.975))
6.5 

ce qui donnera un intervalle d’acceptation à 95%. Dans le second cas, on prendra une région de rejet qui commencera à la note minimal (-10), et on montera jusqu’au quantile à 95% de G, soit
> (note.sup=G.inv(.95))

Moralité, il semblerait que pour 20% des élèves, on accepte l’hypothèse comme quoi ils ont répondu au hasard, et on atteint 30% si on prend l’autre hypothèse alternative. En fait, étrangement, si on regarde les réponses à la première question sur le calcul de la moyenne, presque qu’un tiers des élèves a répondu que s’ils cochaient au hasard, ils auraient 6,6667/20. Or ce tiers a précisément obtenu 6,83/20. Donc il semble qu’un tiers des personnes aient effectivement répondu hasard !
Si j’autorise la personne à choisir (au hasard) l’option de ne pas répondre, je me retrouve avec une loi dite multinomiale. On prendra par la suite la proba de ne pas répondre à la question comme paramètre, afin de voir comme le test évolue avec ce paramètre.
Pour décrire la loi mutinomiale

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/eqn4809-13.png

Bref, là aussi on peut obtenir la fonction de répartition de cette loi, ainsi que sa densité… au moins numériquement, car la somme de lois multinomiales indé
pendantes, ça peut vite devenir un peu compliqué.. Bon, je vais faire du numérique par la suite…
On va commencer par supposer que dans 10% des cas, on ne sait pas répondre, ou plutôt on ne répond rien… Dans ce cas, on note que cela ne change pas trop (au moins sur la région de rejet du test).

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.rejet-rejet-H0-notes-qcm-x-50_m.jpg

Si par contre on suppose que dans 1/4 des cas, on ne répond pas, la région rétrécit,

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.rejet-rejet-H0-notes-qcm-x-50_m.jpg

et encore davantage si on suppose qu’on ne répond pas dans 50% des cas…

https://blogperso.univ-rennes1.fr/arthur.charpentier/public/perso2/.rejet-rejet-H0-notes-qcm-x-50_m.jpg

1 je ne prends pas en compte la possibilité qu’une personne puisse ne pas répondre…. on verra ça plus tard.

Intervalle de confiance pour une proportion

Je voulais poster un petit billet puisque c’est la deuxième fois que l’on me pose une question similaire. Pour déterminer un intervalle de confiance pour une proportion, on connaît tous la formule

La question qui m’avait été posée est simple: que se passe-t-il si n est petit, et si – empiriquement – on a trouvé une proportion de 100% ? Dit autrement: “on fait 5 lancers de “pile” ou “face” et que on obtient 5 fois “pile”: que peut-on dire sur p, la probabilité de tomber sur “pile” ?”

On utilisant la formule précédente, on est un peu bloqué…

  • La première réponse pourrait être fréquentiste, et consiste à proposer une lecture duale de l’intervalle de confiance, et des abaques.

Le graphique de droite montre, en trait violet, l’intervalle de confiance approché, c’est à dire avec l’intervalle de confiance indiqué ci-dessus. La courbe bleue, encadrant la zone verte correspond à l’intervalle de confiance calculé à partir des quantiles d’une loi binomiale de paramètre la fréquence empirique observée. Ici, on aurait envie de dire que p a 95% de chance d’être supérieure à 47.8%.

En fait, cette “lecture duale” des abaques revient à construire – me semble-t-il – l’intervalle de confiance de Clopper-Pearson. L’intervalle de confiance a ici la forme suivante

Notons que c’est cette technique qui est utilisée sous R, via la fonction binom.test du package de base. Ici, on obtient

> binom.test(x, n, p = 0.5,conf.level = 0.95)
Exact binomial test
data: x and n
number of successes = 5, number of trials = 5, p-value = 0.0625
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.4781762 1.0000000
sample estimates:
probability of success
1

…. ce qui correspond exactement à l’intervalle de confiance calculé précédemment. On retrouve d’ailleurs cette fonction dans binconf i.e.

> binconf(x = 5, n = 5, alpha = .05, method = "exact")
PointEst Lower Upper
1 0.4781762 1

qui donne (c’est rassurant) toujours la même chose. On peut aussi programme l’intervalle de confiance proposé par Heldge Blacker, qui donne ici un intervalle de confiance de la forme [50%; 100%]. Pour aller plus loin, je pourrais renvoyer au papier de Lawrence Brown, Tony Cai et Anirban DasGupta, dans Statistical Science sur ce sujet.

  • Une seconde réponse pourrait être une approche bayésienne, qui d’ailleurs correspond exactement à l’expérience de Bayes (1763) ou Laplace (1786).

On suppose que p est une variable aléatoire. On se donne une information a priori, et on peut se demander ce que devient la loi a posteriori de p, sachant que 5 “piles” ont été observés. On suppose que les tirages sont indépendants, et donc la probabilité d’avoir k fois “piles” sur n tirages suit une loi binomiale, conditionnellement à p, i.e.

Le plus simple est de considérer comme loi a priori la loi conjuguée de la loi binomiale, à savoir une loi beta. Dans ce cas, si p suit a priori une loi B(a,b), alors la loi a posteriori de p, sachant que x “pile” ont été observés est une loi B(a+x,b+n-x).

expérience décrite

La courbe en rouge est la densité de la loi a priori (uniforme ou beta de paramètres 2 et 2), et la courbe bleue, la loi a posteriori, avec un intervalle de confiance à 95%. Sur l’exemple de droite, par exemple, on a 95% de chance que p soit compris entre 47,3% et 96,8%.