Tag Archives: ajustement

Test d’ajustement et lâcher de bombes

Avant les applications demain en cours, un petit billet (presque d’actualité) sur une application évoquée vendredi dernier sur le test du chi-deux: l’ajustement de lois.

Le problème est le suivant: un pays se fait bombarder. Et les dirigeants doivent se demander si certaines cibles sont visées (auquel cas il peut être légitime de les déplacer) ou au contraire si les tirs sont aléatoires. Quand je dis que des cibles sont visées, j’entends aussi par là que les pilotes savent viser. Car les pilotes ont (je pense, ou j’espère) un carnet de route à suivre, avec des cibles à viser

Ce que l’on va se demander c’est plutôt si les pilotes savent – ou arrive à – viser. Bref, le problème peut se modéliser en supposant que le lancer de bombes se fait (globalement) selon un processus de Poisson. Localement, si les tirs sont aléatoires, on devrait observer des tirages de lois de Poisson. C’est tout du moins la théorie que R. D. Clarke (alors actuaire chez Prudential Assurance Company) a utilisé pendant la seconde guerre mondiale, sur les bombardements à Londres (en ligne ici). Cette histoire (vraie) est reprise par Thomas Pynchon dans Gravity’s Rainbow (en ligne ici)

(l’exemple est même repris dans Feller). Bref, le problème n’est pas de savoir quel serait le paramètre de la loi de Poisson, mais si la loi de Poisson est adaptée, ou pas. C’est ce qu’on appelle un problème d’ajustement de loi.

  • Test du chi-deux et ajustement de lois

Jusqu’à présent, on avait supposé que les observations suivaient une certaine loi, e.g. une loi de Poisson https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-16.gif, et on cherchait à tester une hypothèse de la forme

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-13.gif

versus

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-14.gif

Ici on va cherche à utiliser un test sur une loi multinomiale, de la forme

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-01.gif

versus

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-02.gif

L’hypothèse nulle est ici une égalité vectorielle,

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-03.gif

ou encore

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-56.gif

Dans le cas d’un test d’ajustement de lois de Poisson, si on suppose que les observations suivent une loi https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-09.gif, on va utiliser un test sur une loi multinomiale, avec

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-08.gif

Le soucis est que, puisque l’on se limite à un vecteur de taille finie, le vecteur ne sera pas dans le simplexe (cf ici). Donc classiquement, pour la dernière valeur, on retient une hypothèse de la forme

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-10.gif

Le test est alors basé sur la statistique de Pearson

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-11.gif

qui va suivre (si effectivement les observations suivent une loi de Poisson https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-09.gif) une loi https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-12.gif.
Rappelons que cette statistique peut aussi s’écrire

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-21.gif

 

  • Application pour étudier la précision d’un lancer de bombes

Pendant la second guerre mondiale, R.D. Clarke étudia les impacts de bombes V1 et V2 tombées dans une zone de 144 km2 dans le sud de Londres (l’article original, publié après guerre en 1946, est en ligne ici). Il divisa cette zone en 576 zones de 0,25 km2 et compta le nombre d’impact dans chacune des zones. Il obtint plus précisément la distribution suivante

nbre impacts par zone 0 1 2 3 4 5 et plus
fréquence (nbre zones) 229 211 93 35 7 1

(en fait, on sait que le “5
et plus” correspond à 7, car sait qu’il y a eu 537 bombes sur 576 zones). Avant de se lancer tête baissée, réfléchissons un peu au type de loi que l’on pourrait utiliser. Pour cela, notons https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/Nb.gif le nombre de points tombés dans un ensemble https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/asubb.gif (ou https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/calA.gif désigne la région globale). Si on suppose qu’un nombre aléatoire https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/Npois.gif de points sont lancés aléatoirement dans https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/calA.gif, et que https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/Npois.gif suit une loi de Poisson https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/lambda.gif, alors https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/Nb.gif suit une loi de Poisson de paramètre

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/ratio-aire.gif

Bref, si on regarde plusieurs plusieurs régions https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/calB.gif (de même taille, éventuellement – afin de garder toutes les observations – formant une partition de https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/calA.gif) et si on observe une loi de Poisson, c’est que dans la région https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/calA.gif, les tirs sont faits au hasard.

> (n=c(229,211,93,35,7,0,0,1))
[1] 229 211  93  35   7   0   0   1
> y=0:7
> (lambda=sum(n*y)/sum(n))
[1] 0.9322917
> prob=dpois(y,lambda)
> freq.theo=sum(n)*prob
> freq.emp =n
> cbind(y,trunc(freq.theo),freq.emp)
     y     freq.emp
[1,] 0 226      229
[2,] 1 211      211
[3,] 2  98       93
[4,] 3  30       35
[5,] 4   7        7
[6,] 5   1        0
[7,] 6   0        0
[8,] 7   0        1

On peut commencer par regarder la log-vraisemblance

> logvrais=function(L){sum(log(dpois(y,L))*n)}
> param=seq(.5,1.5,by=.025)
> LV=sapply(param,logvrais)
> plot(param,LV,type="b",col="blue")
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/logvrais1-bombes.png

La statistique du chi-deux, elle, ressemble à ça

> chi2=function(L)sum(n)*sum(((n/sum(n)-dpois(y,L))^2)/dpois(y,L))
> C2=sapply(param,chi2)
> plot(param,C2,type="b",col="red")xxx.
https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-bombes.png

Bon, le soucis c’est qu’en tronquant le vecteur (on suppose que le nombre maximum d’impact est 7), la somme des probabilités ne fait pas exactement un. On peut regrouper dans une même classe les fréquences élevées (comme le fait Clarke dans le papier initial d’ailleurs).

> (n=c(229,211,93,35,7,1))
[1] 229 211  93  35   7   1
> y=0:4
> prob=c(dpois(y,lambda),1-ppois(4,lambda))
> freq.theo=sum(n)*prob
> freq.emp =n
> (Q=sum(1
[1] 1.169155
> 1-pchisq(Q,length(y)-1)
[1] 0.8831505

On retrouve les quantités évoquées dans l’article de Clarke. La  statistique de test vaut 1.16 et la p-value associée est de l’ordre de 88%. On va donc accepter l’hypothèse de loi de Poisson,

> chi2=function(L){
+ sum(n)*sum(2)^2)/
+ c(dpois(y,L),1-ppois(4,L)) )}
> C2=sapply(param,chi2)
> plot(param,C2,type="b",col="red")

La valeur de la statistique du chi-deux en fonction du paramètre de la loi de Poisson est représentée ci-dessous,

https://f.hypotheses.org/wp-content/blogs.dir/253/files/2015/12/chi2-bombes-v2.png

et le trait horizontal est la valeur seuil de la région critique (en dessous, on accepte l’ajustement d’une loi de Poisson),

> abline(h=qchisq(.95,length(y)-1),lty=2)

Il existe plusieurs fonctions sous R permettant de faire des choses semblables,

> library(vcd)
> (n=c(229,211,93,35,7,0,0,1))
[1] 229 211  93  35   7   0   0   1
> nsim=c(rep(y[0],n[0]),rep(y[1],n[1]),
       rep(y[2],n[2]),rep(y[3],n[3]),
       rep(y[4],n[4]),rep(y[5],n[5]),
       rep(y[6],n[6]),rep(y[7],n[7]),
       rep(y[8],n[8]))
> gf=goodfit(nsim,type="poisson",method="ML")
> summary(gf)
 
 Goodness-of-fit test for poisson distribution
 
 X^2 df P(> X^2)
Likelihood Ratio 4.007784 3 0.2606249
 
> gf=goodfit(nsim,type="poisson",method="MinChisq")
> summary(gf)
 
 Goodness-of-fit test for poisson distribution
 
 X^2 df P(> X^2)
Pearson 1.275499 3 0.7349592

Bref, quelle que soit la méthode utilisée, on notera que l’on accepte toujours l’hypothèse d’une loi de Poisson. Autrement dit les bombes étaient envoyés un peu au hasard dans Londres….

Et si on y réfléchit un peu, d’un point de vue de la théorie des jeux, c’est effectivement une stratégie optimale…

  1. freq.theo-freq.emp)^2)/freq.theo []
  2. n/sum(n)-c(dpois(y,L),1-ppois(4,L []

Statistiques, de la théorie à la pratique

Afin d’illustrer un peu le cours, je mets en ligne ici. Il s’agit de coûts de sinistres ayant dépassé 1,000$ (qui ont été ordonnés de manière croissante). Pour importer la base en R, le code est simplement,

> X=read.table(
"http://freakonometrics.blog.free.fr/public/data/sinistres.txt")$x
> X
 [1] 1.003 1.016 1.023 1.027 1.037 1.039 1.0
[13] 1.061 1.072 1.078 1.082 1.087 1.094 1.0
[25] 1.110 1.112 1.117 1.132 1.138 1.141 1.1
[37] 1.180 1.186 1.187 1.190 1.193 1.203 1.2
[49] 1.321 1.326 1.338 1.342 1.343 1.344 1.3
[61] 1.428 1.432 1.442 1.457 1.463 1.466 1.4
[73] 1.551 1.553 1.566 1.584 1.632 1.695 1.6
[85] 1.881 1.893 1.897 1.958 2.000 2.175 2.2
[97] 3.045 3.103 4.495 5.614

Pour information, l’histogramme ressemble à ça

On va supposer que l’on peut ajuster une loi de Pareto, dont la fonction de survie est donnée par

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-1.gif

et de densité

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-2.gif

L’espérance est

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-3.gif

et la variance

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-5.gif
  • Estimateur de la méthode des moments
http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-6.gif

Numériquement, on obtient

> mean(X)/(mean(X)-1)
[1] 2.946827
  • Estimateur du maximum de vraisemblance

On peut écrire la vraisemblance, et on recherche – numériquement – la maximum (ou le minimum de l’opposé de la log-vraisemblance car la fonction recherche ici les minimum),

> f=function(x,a){
+ a/x^(a+1)
+ }
> LogV = function(a,echantillon){
+ -sum(log(f(echantillon,a)))
+ }
> optim(fn=LogV,par=2,echantillon=X,method="BFGS")$par
[1] 2.845093

Mais on aurait pu faire mieux… car ici le maximum de vraisemblance peut s’écrire analytiquement

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-7.gif
> length(X)/(sum(log(X)))
[1] 2.845093

(tout va bien, on trouve pareil).

  • Estimateur par régression linéaire

On peut noter que

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-8.gif

ce qui veut dire qu’on peut tracer le nuage de points du log des couts, et s’ils sont alignés suivant une droite (passant par l’origine), la pente nous donnera un estimateur de notre paramètre. Le plus simple est de prendre la pente de la droite qui passe par l’origine, et par le centre de gravité du nuage,

> Z=X
> Y=log((length(Z):1)/length(Z))
> X=log(Z)
> plot(X,Y)
> -mean(Y)/mean(X)
[1] 2.753414
  • Estimation à l’aide de la médiane

Enfin, notons que la médiane est

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-9.gif

Un estimateur naturel est alors

http://freakonometrics.blog.free.fr/public/perso2/pareto-2700-11.gif

Numériquement, c’est alors tout simple,

> log(2)/log(quantile(X,.5))
     50% 
2.413719

Moralité, on a quatre estimateurs, qui nous donne quatre valeurs numériques différentes

2.9468 2.8451 2.7534 2.4137

Pour ceux qui veulent trouver d’autres estimateurs (car on peut en construire des centaines), une référence intéressante est la suivante,
Quandt, R. E. (1966). Old and new methods of estimation and the Pareto distribution, Metrika, 10, 55-82 (ici)
Le but de la statistique mathématique est de mieux connaître les propriétés (fondamentales) de ces estimateurs afin de choisir le “meilleur“. Et c’est un problème important, par exemple si on souhaite mettre en place un contrat de réassurance (et transférer les très grands sinistres). Par exemple si on veut transférer les sinistres dépassant 100,000$, la probabilité de dépasser ce seuil est 10 fois plus grande avec le dernier estimateur par rapport au premier…

> 1/100^2.4137
[1] 1.48799e-05
 
> 1/100^2.9468
[1] 1.277615e-06

Moralité, les résultats théoriques que l’on aborde en ce moment vont nous aider à résoudre des questions très pratiques…. à suivre.