Lâcher de bombes et facteur d’échelle

Pour ceux qui auraient pris mon précédant billet (en ligne ici) un peu trop au sérieux (oui, il y en a) deux petites précisions,

  • tout d’abord ce n’est pas bien de lancer des bombes sur des gens que l’on ne connait pas (oui, j’ai eu un courriel sur ce point),
  • ensuite, l’étude que je mentionnais date un peu: il y a 50 ans, les outils statistiques étaient limités….

En fait, on s’en doute, la notion d’échelle est fondamentale. Vu de loin, les tirs ne sont pas du tout aléatoires ! ils sont même relativement précis.

 

 

C’est lorsque l’on regarde de manière beaucoup plus fine que le hasard apparait.

Les données précédentes étaient tirés d’une carte google dont j’ai extrait les latitudes et les longitudes (ici), et la base est malheureusement trop incomplète pour être utile. Mais il existe d’autres sources. Par exemple Charles Franklin a mis en ligne une étude intéressante (), et il a bien voulu m’envoyer les données (qu’il a saisi auparavant à la main). De la même manière que dans le billet précédant, on compte les points dans une grille relativement fine (ici 1km de long et de large) centrée sur Londres. Les données brutes sont à gauche, et les données lissées à droite,

 

 

L’hypothèse d’observations i.i.d. utilisée implicitement dans la précédente étude (mais que l’on ne pouvait pas discuter faute de données) ne semble pas vraiment tenir la route.  Et on voit que les régions ayant eu la plus grosse concentration de bombe sont très proches, ce qui traduit l’idée que les tirs ne sont pas complètement aléatoire. Plus on regarde sur une échelle fine, plus on a cette impression (à l’échelle de la rue, la probabilité d’être touché par une bombe est la même pour tous les bâtiments; mais pas à l’échelle de la région).
Mais pour aller plus loin, il faudrait des notions plus poussées que celles abordées en cours, et encore une fois, le but était juste de proposer une application (réelle) de test d’ajustement de loi par un test du chi-deux.

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 http://freakonometrics.hypotheses.org/files/2015/12/chi2-16.gif, et on cherchait à tester une hypothèse de la forme

http://freakonometrics.hypotheses.org/files/2015/12/chi2-13.gif

versus

http://freakonometrics.hypotheses.org/files/2015/12/chi2-14.gif

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

http://freakonometrics.hypotheses.org/files/2015/12/chi2-01.gif

versus

http://freakonometrics.hypotheses.org/files/2015/12/chi2-02.gif

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

http://freakonometrics.hypotheses.org/files/2015/12/chi2-03.gif

ou encore

http://freakonometrics.hypotheses.org/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 http://freakonometrics.hypotheses.org/files/2015/12/chi2-09.gif, on va utiliser un test sur une loi multinomiale, avec

http://freakonometrics.hypotheses.org/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

http://freakonometrics.hypotheses.org/files/2015/12/chi2-10.gif

Le test est alors basé sur la statistique de Pearson

http://freakonometrics.hypotheses.org/files/2015/12/chi2-11.gif

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

http://freakonometrics.hypotheses.org/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 http://freakonometrics.hypotheses.org/files/2015/12/Nb.gif le nombre de points tombés dans un ensemble http://freakonometrics.hypotheses.org/files/2015/12/asubb.gif (ou http://freakonometrics.hypotheses.org/files/2015/12/calA.gif désigne la région globale). Si on suppose qu’un nombre aléatoire http://freakonometrics.hypotheses.org/files/2015/12/Npois.gif de points sont lancés aléatoirement dans http://freakonometrics.hypotheses.org/files/2015/12/calA.gif, et que http://freakonometrics.hypotheses.org/files/2015/12/Npois.gif suit une loi de Poisson http://freakonometrics.hypotheses.org/files/2015/12/lambda.gif, alors http://freakonometrics.hypotheses.org/files/2015/12/Nb.gif suit une loi de Poisson de paramètre

http://freakonometrics.hypotheses.org/files/2015/12/ratio-aire.gif

Bref, si on regarde plusieurs plusieurs régions http://freakonometrics.hypotheses.org/files/2015/12/calB.gif (de même taille, éventuellement – afin de garder toutes les observations – formant une partition de http://freakonometrics.hypotheses.org/files/2015/12/calA.gif) et si on observe une loi de Poisson, c’est que dans la région http://freakonometrics.hypotheses.org/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")
http://freakonometrics.hypotheses.org/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.
http://freakonometrics.hypotheses.org/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,

http://freakonometrics.hypotheses.org/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 []

“Je comprends pas tout, mais j’aime bien”

http://freakonometrics.blog.free.fr/public/perso2/.2249021524.4_s.jpgPlusieurs fois, j’ai lu des commentaires sur le blog (ou sur le premier, à Rennes 1) commençant par quelque chose du genre “j’ai bien aimé même si j’ai pas tout compris“…
Étrangement, je pense à cette phrase régulièrement depuis quelques mois. En fait, en novembre dernier, j’avais acheté pour mon fils un livre écrit par Ivar Ekeland, Le chat au pays des nombres. S’il a bien aimé, c’est sans commune mesure avec la passion que voue ma fille pour ce livre ! Je dois lui lire au moins une fois par semaine ! Elle ne s’en lasse pas ! Et moi non plus d’ailleurs (je finirais par mieux connaitre ce livre que celui de théorie des jeux qu’il avait publié il y a un peu plus longtemps).

http://freakonometrics.blog.free.fr/public/perso2/.1296324_3197388_s.jpg http://freakonometrics.blog.free.fr/public/perso2/.1296324_3197389_s.jpg
Pourtant, le livre explique tout simplement que http://freakonometrics.blog.free.fr/public/perso2/ensembleQ.gif est dénombrable, en construisant une bijection entre http://freakonometrics.blog.free.fr/public/perso2/ensembleN.gif et http://freakonometrics.blog.free.fr/public/perso2/ensembleNN.gif, comme l’avait Cantor (et que l’infini, c’est quand même super grand).
http://freakonometrics.blog.free.fr/public/perso2/cantor_set.jpg

Et comme c’est un livre pour enfants, Ivar utilise l’approche proposée par David Hilbert (ici), à savoir l’hôtel infini. Chaque fraction (i.e. un élément de http://freakonometrics.blog.free.fr/public/perso2/ensembleQ.gif, en bleu au dessus) a une place dans une chambre de l’hôtel infini (dont le numéro de chambre est dans http://freakonometrics.blog.free.fr/public/perso2/ensembleN.gif, en rouge au dessus).
Le livre est génial, et je pousse tout le monde à le livre ! Mais je ne cesse de me demander ce que ma fille en retient, ce qu’elle comprend vraiment…  Car Ivar Ekeland est formidable: il arrive à parler de concept non triviaux à des enfants de 5 ans (pour les plus grands, je peux renvoyer ici pour une relecture du papier de John Nash sur l’interprétation d’une théorème de point fixe).
(Je pourrais aussi noter que ce livre me fait aussi m’interroger sur la notion d’impact factor des publications des chercheurs. Ivar est connu dans la communauté scientifique par bon nombre de papiers théoriques publiés dans les revues les plus prestigieuses, mais je doute qu’un chercheur ait passé autant de temps sur ses articles que ma fille sur son livre.)

Cochrane, Pearson et le test du chi-deux

En cours, nous avons poursuivi sur la loi multinomiale, et le test du chi-deux. Je voulais mettre un petit billet pour récapituler les différents points, et montrer une application numérique (nous reviendrons en détails mercredi sur des applications des outils vus jusqu’à présent).

  • Inférence avec la loi multinomiale

On suppose qu’une variable http://freakonometrics.blog.free.fr/public/maths/coch-01.gif peut prendre http://freakonometrics.blog.free.fr/public/maths/coch-02.gif modalités, notées http://freakonometrics.blog.free.fr/public/maths/coch-03.gif, avec http://freakonometrics.blog.free.fr/public/maths/coch-04.gif. On posera

http://freakonometrics.blog.free.fr/public/maths/coch-05.gif

en notant que http://freakonometrics.blog.free.fr/public/maths/coch-06.gif appartient au simplexe de http://freakonometrics.blog.free.fr/public/maths/coch-07.gif au sens où

http://freakonometrics.blog.free.fr/public/maths/coch-08.gif

On a vu que l’estimateur du maximum de vraisemblance s’obtenait en faisant un peu d’optimisation sous contrainte, et que

http://freakonometrics.blog.free.fr/public/maths/coch-10.gif

(en reprenant les notations du cours). On avait montré que

http://freakonometrics.blog.free.fr/public/maths/coch-11.gif

et on a vu

http://freakonometrics.blog.free.fr/public/maths/coch-13.gif

(ce qui peut se retrouver en introduisant la variable binomiale http://freakonometrics.blog.free.fr/public/maths/coch-16.gif). Mais plus généralement, on finira les calculs permettant d’établir que, pour http://freakonometrics.blog.free.fr/public/maths/coch-17.gif

http://freakonometrics.blog.free.fr/public/maths/coch-18.gif

(ce qui permet d’obtenir la matrice de variance covariance de http://freakonometrics.blog.free.fr/public/maths/coch-20.gif).
En utilisant le théorème central limite on peut montrer que

http://freakonometrics.blog.free.fr/public/maths/coch-23.gif

Sous une forme multivariée, on écrira

http://freakonometrics.blog.free.fr/public/maths/coch-25.gif
http://freakonometrics.blog.free.fr/public/maths/coch-26.gif avec http://freakonometrics.blog.free.fr/public/maths/coch-27.gif et pour http://freakonometrics.blog.free.fr/public/maths/coch-17.gif, http://freakonometrics.blog.free.fr/public/maths/coch-28.gif.

On notera que la somme de la ième colonne de http://freakonometrics.blog.free.fr/public/maths/coch-29.gif est http://freakonometrics.blog.free.fr/public/maths/coch-30.gif.
Aussi, http://freakonometrics.blog.free.fr/public/maths/coch-29.gif n’est pas inversible (c’est le fait que notre paramètre appartient au simplexe).
Pour s’en sortir, la première idée est d’utiliser un peu d’algèbre linéaire. Une matrice http://freakonometrics.blog.free.fr/public/maths/coch-31.gif est une matrice de projection si elle est idempotente, i.e. http://freakonometrics.blog.free.fr/public/maths/coch-32.gif. Ses valeurs propres sont alors 0 ou 1, et si http://freakonometrics.blog.free.fr/public/maths/coch-35.gif est le nombre de fois où 1 est valeur propre, et si http://freakonometrics.blog.free.fr/public/maths/coch-36.gif, alors http://freakonometrics.blog.free.fr/public/maths/coch-37.gif.
Posons http://freakonometrics.blog.free.fr/public/maths/coch-38.gif. Alors

http://freakonometrics.blog.free.fr/public/maths/coch-39.gif

Or compte tenu de la forme de http://freakonometrics.blog.free.fr/public/maths/coch-29.gif,

http://freakonometrics.blog.free.fr/public/maths/coch-40.gif

qui est une matrice de projection dont la trace est http://freakonometrics.blog.free.fr/public/maths/coch-41.gif (qui est aussi le nombre de fois où 1 est valeur propre). Donc

http://freakonometrics.blog.free.fr/public/maths/coch-42.gif

Le test du chi-deux sera basé sur le fait qu’asymmptotiquement

http://freakonometrics.blog.free.fr/public/maths/coch-44.gif

Une autre idée consiste à construire http://freakonometrics.blog.free.fr/public/maths/coch-41.gif variables aléatoires qui seront indépendantes. Mais on peut plutôt regarder les applications de ce test, en particulier comme test d’indépendance.
Pour information, Frank Yates a proposé un correction “pour continuité“, ici. La statistique considérée est alors

http://fre<br /><br /><br /><br />
akonometrics.blog.free.fr/public/maths/coch-46.gif
  • Retour sur le théorème de Cochrane

Soit http://freakonometrics.blog.free.fr/public/maths/coch-50.gif de dimension http://freakonometrics.blog.free.fr/public/maths/coch-51.gif. Posons http://freakonometrics.blog.free.fr/public/maths/coch-59.gif, pour http://freakonometrics.blog.free.fr/public/maths/coch-04.gif, où on notera http://freakonometrics.blog.free.fr/public/maths/coch-60.gif le rang de http://freakonometrics.blog.free.fr/public/maths/coch-62.gif, en supposant que les http://freakonometrics.blog.free.fr/public/maths/coch-62.gif sont positives semidefinies, alors on a équivalence entre

  •  http://freakonometrics.blog.free.fr/public/maths/coch-63.gif
  • http://freakonometrics.blog.free.fr/public/maths/coch-64.gif pour http://freakonometrics.blog.free.fr/public/maths/coch-04.gif,
  • les http://freakonometrics.blog.free.fr/public/maths/coch-65.gif sont des variables indépendantes.

Les http://freakonometrics.blog.free.fr/public/maths/coch-65.gif s’interprètent comme des longueurs (euclidienne) de projections d’un vecteur Gaussien sur des sous-espaces orthogonaux (de dimension respective http://freakonometrics.blog.free.fr/public/maths/coch-60.gif). Si ces vecteurs sont indépendants, et suivent des lois du chi-deux à http://freakonometrics.blog.free.fr/public/maths/coch-60.gif degrés de libertés, avec http://freakonometrics.blog.free.fr/public/maths/coch-63.gif, alors les sous-espaces sont orthogonaux, et supplémentaires. On peut y voir une espèce d’extension du théorème de Pythagore, en remplaçant la notion de vecteurs orthogonaux par des variables indépendantes suivant des lois du chi-deux, et la somme des carrés des longueurs devient la somme des degrés de liberté.

  • Application comme test d’indépendance

Considérons deux variables http://freakonometrics.blog.free.fr/public/maths/coch-66.gif pouvant prendre toutes les deux deux modalités (disons deux lois binomiales). On est alors face a une loi multinomiale à 4 modalités

  • http://freakonometrics.blog.free.fr/public/maths/coch-79.gif avec probabilité http://freakonometrics.blog.free.fr/public/maths/coch-70.gif
  • http://freakonometrics.blog.free.fr/public/maths/coch-78.gif avec probabilité http://freakonometrics.blog.free.fr/public/maths/coch-73.gif
  • http://freakonometrics.blog.free.fr/public/maths/10gif.gif avec probabilité http://freakonometrics.blog.free.fr/public/maths/coch-74.gif
  • http://freakonometrics.blog.free.fr/public/maths/coch-77.gif avec probabilité http://freakonometrics.blog.free.fr/public/maths/coch-75.gif

Un test d’indépendance revient à tester si la loi multinomiale peut s’écrire

http://freakonometrics.blog.free.fr/public/maths/chi2-ab.gif
http://freakonometrics.blog.free.fr/public/maths/chi_ab2.gif
http://freakonometrics.blog.free.fr/public/maths/chi_ab3.gif
http://freakonometrics.blog.free.fr/public/maths/chi_ab4.gif

pour des vecteurs http://freakonometrics.blog.free.fr/public/maths/chi-a.gif et http://freakonometrics.blog.free.fr/public/maths/chi-b.gif. tels que http://freakonometrics.blog.free.fr/public/maths/chi-a2.gif et http://freakonometrics.blog.free.fr/public/maths/chi-b2.gif. On a alors http://freakonometrics.blog.free.fr/public/maths/chi212121.gif contraintes sur les paramètres. Ces deux contraintes additionnelles font que la statistique de test s’écrit

http://freakonometrics.blog.free.fr/public/maths/CHI-INDEP.gif

qui va suivre asymptotiquement une loi http://freakonometrics.blog.free.fr/public/maths/CHI1.gif i.e. http://freakonometrics.blog.free.fr/public/maths/CHI12.gif d’après le théorème de Cochrane.

  • Peine de mort pour les condamnées pour meurtre en Floride 1976-1987

en fonction de la “race” du meurtrier et de la victime,

  • meurtrier de “race blanche” et victime de “race blanche“: 53 condamnés à mort, 414 non condamnés à mort
  • meurtrier de “race blanche” et victime de “race noire“: 0 condamné à mort, 16 non condamnés à mort
  • meurtrier de “race noire“et victime de “race blanche“: 11 condamnés à mort, 37 non condamnés à mort
  • meurtrier de “race noire“et victime de “race noire“: 4 condamnés à mort, 139 non condamnés à mort

On peut alors faire des tests d’indépendance, entre la “race” du meurtrier et le verdict par exemple.

MEURTRIER=matrix(c(53+0,11+4,414+16,139+37),2,2)
VICTIME  =matrix(c(53+11,0+4,414+37,139+16),2,2)
n=sum(MEURTRIER)
(PROBMEUTR=MEURTRIER/n)
           [,1]      [,2]
[1,] 0.07863501 0.6379822
[2,] 0.02225519 0.2611276

SL=rowSums(PROBMEUTR)
SC=colSums(PROBMEUTR)
(MEUTRINDEP=outer(SL, SC, "*"))
           [,1]      [,2]
[1,] 0.07229966 0.6443176
[2,] 0.02859055 0.2547922

(Q=n*sum((PROBMEUTR - MEUTRINDEP)^2/MEUTRINDEP))
[1] 1.468519

(Qcorrec=n*sum((abs(PROBMEUTR - MEUTRINDEP)-.5/n)^2/MEUTRINDEP))
[1] 1.144741

pchisq(Qcorrec, (2-1)*(2-1), lower.tail
 = FALSE)
[1] 0.2846528

qchisq(.95, (2-1)*(2-1))
[1] 3.841459

chisq.test(MEURTRIER)

Pearson's Chi-squared test with Yates' continuity correction

data:  MEURTRIER 
X-squared = 1.1447, df = 1, p-value = 0.2847

On rejette donc l’hypothèse d’indépendance.

911, jour après jour

Après deux billets (ici puis ) sur les cycles intrajournaliers des appels au 911, on peut se demander comment les crimes évoluent au cours de la semaine.

Pour l’ensemble des appels passés au 911, on a la distribution suivante

i.e. un pic les vendredi soir et samedi soir, et un creux le dimanche. Si on regarde les appels pour des cambriolages, on a la distribution suivante

avec des pics en matinée, les vendredi après midi, et les fins de semaine. On peut aussi suivre les troubles de la paix,

qui surviennent certes vers minuit, mais essentiellement en fin de semaine. Ce qui contraste assez avec les hold-ups,

Manifestement, il y a des tendances assez claires. La prochaine étape sera de regarder un peu les saisons, ou mieux, l’impact du climat…

à suivre donc…

Circular or spherical data, and density estimation

I few years ago, while I was working on kernel based density estimation on compact support distribution (like copulas) I went through a series of papers on circular distributions. By that time, I thought it was something for mathematicians working on weird spaces…. but during the past weeks, I saw several potential applications of those estimators.

  • circular data density estimation

Consider the density of an angle say, i.e. a function http://freakonometrics.hypotheses.org/files/2015/12/circ-01.gif such that

http://freakonometrics.hypotheses.org/files/2015/12/circ-02.gif

with a circular relationship, i.e. http://freakonometrics.hypotheses.org/files/2015/12/circ-03.gif. It can be seen as an invariance by rotation.
von Mises proposed a parametric model in 1918 (see here or there), assuming that

http://freakonometrics.hypotheses.org/files/2015/12/circ-04.gif

where http://freakonometrics.hypotheses.org/files/2015/12/circ-05.gif is Bessel modified function of order 1,

http://freakonometrics.hypotheses.org/files/2015/12/circ-06.gif

(which is simply a normalization parameter). There are two parameters here, http://freakonometrics.hypotheses.org/files/2015/12/circ-07.gif (some concentration parameter) and mu a direction.
From a series of observed angleshttp://freakonometrics.hypotheses.org/files/2015/12/circ-08.gif, the maximum likelihood estimator for kappa is solution of

http://freakonometrics.hypotheses.org/files/2015/12/circ-09.gif

where

http://freakonometrics.hypotheses.org/files/2015/12/circ-10.gif

and

http://freakonometrics.hypotheses.org/files/2015/12/circ-11.gif

and where http://freakonometrics.hypotheses.org/files/2015/12/circ-12.gif, where those functions are modified Bessel functions. Well, that estimator is biased, but it is possible to improve it (see here or there). This can be done easily in R (actually Jeff Gill – here – used that package in several applications). But I am not a big fan of that technique….

  • density estimation for hours on simulated data

A nice application can be on the estimation of the daily density of a temporal events (e.g. phone calls as we’ll see later on, or email arrival time). Let http://freakonometrics.hypotheses.org/files/2015/12/circ-13.gif is the time (in hours) for the http://freakonometrics.hypotheses.org/files/2015/12/circ-14.gifth observation (the http://freakonometrics.hypotheses.org/files/2015/12/circ-14.gifth phone call received). Then set

http://freakonometrics.hypotheses.org/files/2015/12/circ-15.gif

The time is now seen as an angle. It is possible to consider the equivalent of an histogram,

set.seed(1)
library(circular)
X=rbeta(100,shape1=2,shape2=4)*24
Omega=2*pi*X/24
Omegat=2*pi*trunc(X)/24
H=circular(Omega,type="angle",units="radians",rotation="clock")
Ht=circular(Omegat,type="angle",units="radians",rotation="clock")
plot(Ht, stack=FALSE, shrink=1.3, cex=1.03,
axes=FALSE,tol=0.8,zero=c(rad(90)),bins=24,ylim=c(0,1))
points(Ht, rotation = "clock", zero =c(rad(90)),
col = "1", cex=1.03, stack=TRUE )

rose.diag(Ht-pi/2,bins=24,shrink=0.33,xlim=c(-2,2),ylim=c(-2,2),
axes=FALSE,prop=1.5)

or a kernel based estimation of the density (the gray line on the right).

circ.dens = density(Ht+3*pi/2,bw=20)
plot(Ht, stack=TRUE, shrink=.35, cex=0, sep=0.0,
axes=FALSE,tol=.8,zero=c(0),bins=24,
xlim=c(-2,2),ylim=c(-2,2), ticks=TRUE, tcl=.075)
lines(circ.dens, col="darkgrey", lwd=3)
text(0,0.8,"24", cex=2); text(0,-0.8,"12",cex=2);
text(0.8,0,"6",cex=2); text(-0.8,0,"18",cex=2)

The code looks rather simple. But I am not very comfortable using codes that I do not completely understand. So I did my own. The first step was to get a graph similar to the one we have on the right, except that I prefer my own kernel based estimator. The idea is that instead of estimating the density on http://freakonometrics.hypotheses.org/files/2015/12/Xi.gif, we estimate it on the sample http://freakonometrics.hypotheses.org/files/2015/12/circular-density-3.gif. Then we multiply by 3 to get the density only on http://freakonometrics.hypotheses.org/files/2015/12/0-24.gif. For the bandwidth, I took the same as the one that we would have taken on http://freakonometrics.hypotheses.org/files/2015/12/Xi.gif

The code is simply the following

U=seq(0,1,by=1/250)
O=U*2*pi
U12=seq(0,1,by=1/24)
O12=U12*2*pi
X=rbeta(100,shape1=2,shape2=4)*24
OM=2*pi*X/24
XL=c(X-24,X,X+24)
d=density(X)
d=density(XL,bw=d$bw,n=1500)
I=which((d$x>=6)&(d$x<=30))
Od=d$x[I]/24*2*pi-pi/2
Dd=d$y[I]/max(d$y)+1

plot(cos(O),-sin(O),xlim=c(-2,2),ylim=c(-2,2), type="l",axes=FALSE,xlab="",ylab="") for(i in pi/12*(0:12)){ abline(a=0,b=tan(i),lty=1,col="light yellow")} segments(.9*cos(O12),.9*sin(O12),1.1*cos(O12),1.1*sin(O12)) lines(Dd*cos(Od),-Dd*sin(Od),col="red",lwd=1.5) text(.7,0,"6"); text(-.7,0,"18") text(0,-.7,"12"); text(0,.7,"24") R=1/24/max(d$y)/3+1 lines(R*cos(O),R*sin(O),lty=2)

Note that it is possible to stress more (visually) on hours having few phone calls, or a lot (compared with an homogeneous Poisson process), e.g.

plot(cos(O),-sin(O),xlim=c(-2,2),ylim=c(-2,2),
type="l",axes=FALSE,xlab="",ylab="")
for(i in pi/12*(0:12)){
abline(a=0,b=tan(i),lty=1,col="light yellow")}
segments(2*cos(O12),2*sin(O12),1.1*cos(O12),1.1*sin(O12), col="light grey")
segments(.9*cos(O12),.9*sin(O12),1.1*cos(O12),1.1*sin(O12))
text(.7,0,"6")
text(-.7,0,"18")
text(0,-.7,"12")
text(0,.7,"24")
R=1/24/max(d$y)/3+1
lines(R*cos(O),R*sin(O),lty=2)
AX=R*cos(Od);AY=-R*sin(Od)
BX=Dd*cos(Od);BY=-Dd*sin(Od)
COUL=rep("blue",length(AX))
COUL[R<Dd]="red"
CM=cm.colors(200)
a=trunc(100*Dd/R)
COUL=CM[a]
segments(AX,AY,BX,BY,col=COUL,lwd=2)
lines(Dd*cos(Od),-Dd*sin(Od),lwd=2)

We get here those two graphs,

To be honest, I do not really like that representation – even if it looks nice. If we compare that circular representation to a more classical one (from 0:00 till 23:59 one the graph on the left, below), I do have a problem to interpret the areas in blue and pink.

density of wind direction

On the left, we compare two densities, so the area in pink is the same as the area in blue. But here, it is no longer the case: the area in pink is always larger to the one in blue. So it might help so see when we have a difference, but there is a scaling issue that we cannot discuss further… But less us see if we can use that estimation technique to several problems.

A standard application when studying angles is wind direction. For instance, in Montréal, it is possible to find hourly observations, starting in 1974 (we just need a R robot to pick up the information, but I’ll tell more about that in another post, someday). Here, we have directly an angle. So we can use a code rather similar to the one used above to estimate the distribution of wind direction in Montréal.

density of 911 phone calls

Note that our estimate is consistent with several graphs that can be found on meteorological websites (e.g. the one above on the right, that was found here).

In a recent post (here) I wanted to check about the “midnight crime” myth, using hours of 911 phone calls in Montréal.

That was for all phone calls. But if we look more specifically, for burglaries, we have the distribution on the left, and for conflicts the one on the right

We do clearly observe that gun shots occur a bit before midnight. See also here for another study, but this time in NYC (thanks @PAC for the link).while for gun shots, we have the distribution on the left, and for “troubles” (basically people making too much noisy in parties) or “noise” the one on the right

  • density of earth temperatures, or earthquakes

Of course it is also possible to work in higher dimension. Before, we went from densities on http://freakonometrics.hypotheses.org/files/2015/12/circ-16.gif to densities on the unit circle http://freakonometrics.hypotheses.org/files/2015/12/circ-18.gif. But similarly, it is possible to go from http://freakonometrics.hypotheses.org/files/2015/12/circ-17.gif to the unit sphere http://freakonometrics.hypotheses.org/files/2015/12/circ-19.gif. A nice application being global climate studies,

The idea being that point on the left above are extremely close to the one on the right. An application can be e.g. on earthquakes occurrence. Data can be found here.

library(ks)
X=cbind(EQ$Longitude,EQ$Latitude)
Hpi1 = Hpi(x = X)
DX=kde(x = X, H = Hpi1)
library(maps)
map("world")
plot(DX,add=TRUE,col="red")
points(X,cex=.2,col="blue")
Y=rbind(cbind(X[,1],X[,2]),cbind(X[,1]+360,X[,2]),
cbind(X[,1]-360,X[,2]),cbind(X[,1],X[,2]+180),
cbind(X[,1]+360,X[,2]+180),cbind(X[,1]-360,X[,2]+180), cbind(X[,1],X[,2]-180),cbind(X[,1]+360, X[,2]-180),cbind(X[,1]-360,X[,2]-180)) DY=kde(x = Y, H = Hpi1) library(maps) plot (DY,add=TRUE,col="purple")

Without any correction, we get the red level curves. The pink one integrates correction.

Want to say one thing and the exact oppositive with strong confidence ?

No need to do politics. Just take a statistical course. And I do not talk about misinterpretation of statistics, but I talk about the mathematical foundations of statistical tests.
Consider the following parametric test, with a one-dimensional parameter: http://freakonometrics.blog.free.fr/public/perso2/test-lies-01.gif versus http://freakonometrics.blog.free.fr/public/perso2/test-lies-02.gif, for some (fixed) http://freakonometrics.blog.free.fr/public/perso2/test-lies-03.gif. A standard way of doing such a test is to consider an rejection region http://freakonometrics.blog.free.fr/public/perso2/test-lies-05.gif. The test works as follows: consider a sample http://freakonometrics.blog.free.fr/public/perso2/test-lies-06.gif,

  • if http://freakonometrics.blog.free.fr/public/perso2/test-lies-07.gif, then we accept http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif
  • if http://freakonometrics.blog.free.fr/public/perso2/test-lies-09.gif, the we reject http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif

For instance, consider the case of a Bernoulli sample, with probability http://freakonometrics.blog.free.fr/public/perso2/test-lies-62.gif. The standard idea is to define

http://freakonometrics.blog.free.fr/public/perso2/test-lies-13.gif

The rejection region is then based on statistic http://freakonometrics.blog.free.fr/public/perso2/test-lies-210.gif,

  • if http://freakonometrics.blog.free.fr/public/perso2/test-lies-25.gif, then we accept http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif
  • if http://freakonometrics.blog.free.fr/public/perso2/test-lies-22.gif, the we reject http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif

where threshold http://freakonometrics.blog.free.fr/public/perso2/test-lies-26.gif is taken so that the probability to make a first type error is http://freakonometrics.blog.free.fr/public/perso2/test-lies-28.gif(say 5%) using the Gaussian approximation for z. Here

http://freakonometrics.blog.free.fr/public/perso2/test-lies-30.gif

Thus, the acceptation region is then the green area below, while the rejection region is the red one, for http://freakonometrics.blog.free.fr/public/perso2/test-lies-210.gif.

Consider now the exact opposite test (with the same http://freakonometrics.blog.free.fr/public/perso2/test-lies-03.gif), http://freakonometrics.blog.free.fr/public/perso2/test-lies-51.gifversus http://freakonometrics.blog.free.fr/public/perso2/test-lies-52.gif. Here, we use the same statistics, and the test is

  • if http://freakonometrics.blog.free.fr/public/perso2/test-lies-22.gif, then we accept http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif
  • if http://freakonometrics.blog.free.fr/public/perso2/test-lies-25.gif, the we reject http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif

where now

http://freakonometrics.blog.free.fr/public/perso2/test-lies-50.gif

Thus, now, the acceptation region is then the green area below, while the rejection region is the red one.

So if we summarize what we just said,

  • in the region on the left below, both test agree that http://freakonometrics.blog.free.fr/public/perso2/test-lies-55.gif
  • in the region on the right below, both test agree that http://freakonometrics.blog.free.fr/public/perso2/test-lies-57.gif
  • and in the region in blue, in the middle, the two tests disagree (one claims that http://freakonometrics.blog.free.fr/public/perso2/test-lies-55.gif, and the other one that http://freakonometrics.blog.free.fr/public/perso2/test-lies-57.gif)

Here is the evolution of the region as a function of http://freakonometrics.blog.free.fr/public/perso2/test-lies-56.gif (the size of the sample) when the sample frequency is 20%. With a small sample size, we can hardly say anything.

n=seq(1,100)
p=0.2
x1=p+qnorm(.95)*sqrt(p*(1-p)/n)
x2=p+qnorm(.05)*sqrt(p*(1-p)/n)
plot(n,x1,type="l",ylim=c(0,1))
polygon(c(n,rev(n)),c(x1,rev(x2)),col="light blue",border=NA)
lines(n,x1,lwd=2,col="red")
lines(n,x2,lwd=2,col="red")

One might say that those bounds are based on a Gaussian approximation which is not correct when http://freakonometrics.blog.free.fr/public/perso2/test-lies-56.gif is too small. So we can compute exact bounds,
y1=qbinom(.95,size=n,prob=p)/n
y2=qbinom(.05,size=n,prob=p)/n
polygon(c(n,rev(n)),c(y1,rev(y2)),col="blue",border=NA)
lines(n,y1,lwd=2,col="red")
lines(n,y2,lwd=2,col="red")

and we get

This is what we can observe if we use R statistical procedures, either the asymptotic one,

> prop.test(2,10,.5,alternative="less")
 
1-sample proportions test with continuity correction
 
data:  2 out of 10, null probability 0.5
X-squared = 2.5, df = 1, p-value = 0.05692
alternative hypothesis: true p is less than 0.5
95 percent confidence interval:
0.0000000 0.5100219
sample estimates:
p
0.2
 
> prop.test(2,10,.5,alternative="greater")
 
1-sample proportions test with continuity correction
 
data:  2 out of 10, null probability 0.5
X-squared = 2.5, df = 1, p-value = 0.943
alternative hypothesis: true p is greater than 0.5
95 percent confidence interval:
0.04368507 1.00000000
sample estimates:
p
0.2

or a more accurate one

> binom.test(2,10,.5,alternative="less")
 
Exact binomial test
 
data:  2 and 10
number of successes = 2, number of trials = 10, p-value = 0.05469
alternative hypothesis: true probability of success is less than 0.5
95 percent confidence interval:
0.0000000 0.5069013
sample estimates:
probability of success
0.2
 
> binom.test(2,10,.5,alternative="greater")
 
Exact binomial test
 
data:  2 and 10
number of successes = 2, number of trials = 10, p-value = 0.9893
alternative hypothesis: true probability of success is greater than 0.5
95 percent confidence interval:
0.03677144 1.00000000
sample estimates:
probability of success
0.2

Here, when the sample frequency is 20% and http://freakonometrics.blog.free.fr/public/perso2/test-lies-56.gif is equal to 10, we accept at the same time that theta is higher than 50% and lower than 50%.
And obviously it is not only a theoretical problem: it has obviously some strong implications. This morning, a good friend mentioned a post published some months ago, online here, about discrimination, and the lack of women with academic positions in mathematics, in France. As claimed by the author of the post“A Paris VI, meilleure université française selon son président, sur 11 postes de maitres de conférences, 5 filles classées premières. Il y a donc des filles excellentes ? A Toulouse, sur 4 postes, 2 filles premières. Parité parfaite. Mais à côté de cela, Bordeaux, 4 postes, 0 fille première. Littoral, 3 postes, 0 fille, Nice, 5 postes, 0 fille, Rennes, 7 postes, 0 fille…”.
Consider the latter one: in Rennes, out of 7 people hired last year, no woman. So in some sense, it looks obvious that there is some kind of discrimination ! Zero out of seven ! Well, if we consider the fact that around 30% of PhD thesis in mathematics were defended by women those years, we can also try to see is there if no “positive discrimination“, i.e. test http://freakonometrics.blog.free.fr/public/perso2/test-lies-60.gif where theta is the probability to hire a woman (just to be a little bit provocative).

> prop.test(0,7,.3,alternative="less")
 
1-sample proportions test with continuity correction
 
data:  0 out of 7, null probability 0.3
X-squared = 1.7415, df = 1, p-value = 0.09347
alternative hypothesis: true p is less than 0.3
95 percent confidence interval:
0.0000000 0.3719021
sample estimates:
p
0
 
Warning message:
In prop.test(0, 7, 0.3, alternative = "less") :
Chi-squared approximation may be incorrect
> binom.test(0,7,.3,alternative="less")
 
Exact binomial test
 
data:  0 and 7
number of successes = 0, number of trials = 7, p-value = 0.08235
alternative hypothesis: true probability of success is less than 0.3
95 percent confidence interval:
0.0000000 0.3481637
sample estimates:
probability of success
0

With no woman hired that year, we can still pretend that there was some kind of “positive discrimination“. An note that we do accept – with more confidence – the assumption of “positive discrimination” if we look at all universities together,

> prop.test(5+2,11+4+4+3+5+7,.3,alternative="less")
 
1-sample proportions test with continuity correction
 
data:  5 + 2 out of 11 + 4 + 4 + 3 + 5 + 7, null probability 0.3
X-squared = 1.021, df = 1, p-value = 0.1561
alternative hypothesis: true p is less than 0.3
95 percent confidence interval:
0.0000000 0.3556254
sample estimates:
p
0.2058824
 
> binom.test(5+2,11+4+4+3+5+7,.3,alternative="less")
 
Exact binomial test
 
data:  5 + 2 and 11 + 4 + 4 + 3 + 5 + 7
number of successes = 7, number of trials = 34, p-value = 0.1558
alternative hypothesis: true probability of success is less than 0.3
95 percent confidence interval:
0.0000000 0.3521612
sample estimates:
probability of success
0.2058824

So obviously, with small sample, almost anything can be claimed !

Du sex-ratio en France

En novembre dernier, Baptiste @Coulmont m’avait envoyé un courriel correspondant à ce qu’il a mis en ligne sur son blog (ici) sur l’utilisation du fichier des prénoms français pour analyser le sex ratio à la naissance en France. J’attendais qu’il publie ses commentaires avant de mettre les miens (car le graphique qu’il a mis en ligne ce matin m’avait fait m’interroger).
Pour faire simple, la base de prénoms inclue un sexe. Mais si on regarde le rapport du nombre de garçon sur le nombre de filles, à la naissance, on a le graphique ci-dessous,

dat=read.table("nat2004.csv",sep=";",header=TRUE)
naissancesm=rep(NA,105)
for (i in 1900:2004) {
naissancesm[i-1899]=sum(dat[dat$annais==i&dat$sexe==1,"nombre"],
na.rm=TRUE)
}
naissancesf=rep(NA,105)
for (i in 1900:2004) {
naissancesf[i-1899]=sum(dat[dat$annais==i&dat$sexe==2,"nombre"],
na.rm=TRUE)
}
plot(1900:2004,naissancesm/naissancesf,col="red")

La tendance du début est très surprenante. Par exemple,  si on reprend les chiffres donnés par Pierre Simon Laplace sur les naissances à Paris entre 1750 et 1800 (mentionné hier, ici) on est déjà sur ratio de l’ordre de 1.05 (que l’on retrouve sur la fin de notre graphique, mais pas le début).
> 393386/377555
[1] 1.041930
 
> 251527/241945
[1] 1.039604

Donc il n’y a pas de raison d’avoir cette diff1rence. La conclusion semble être qu’il y a un soucis sur la base des prénoms. En effet, dans un rapport de Anouch Chahnazarian (ici), on retrouve l’évolution suivante,

que l’on peut rapprocher des données d’Éric Brian & Marie Jaisson (ici pour les données et quelques pages) utilisées dans leur ouvrage Le sexisme de la première heure, hasard et sociologie, qui mesure ici la fréquence de garçons à la naissance

Si l’on compare ces dernières données (via le fichier ici), calculé sur les données de l’INSEE et de l’INED, on retrouve un niveau très proche de celui que l’on a sur le fichier des prénoms, avec toutefois un biais constamment négatif.

b=read.table("http://freakonometrics.blog.free.fr/public/data/sex-ratio.txt")
X=b$V1
Y=b$V2/(100-b$V2)
lines(X,Y,col="blue")

Attention donc à la variable de sexe dans la base surtout avant guerre (en espérant que ce soit le seul soucis). Et je renvoie au blog de Baptiste qui publie toujours des choses très amusantes sur les prénoms, ici (je n’ai plus trop eu le temps de travailler dessus depuis les 5 billets en ligne ).

Minuit, l’heure du crime ?

Il y a quelques jours, j’étais invité par Rémi pour me rendre au poste de police pour discuter de la possibilité de consulter quelques données. Et Rémi m’a donné un historique de près d’un million d’appels téléphonique d’urgence, sur Montréal, passés au cours des six dernières années. Histoire de prendre un peu en main les données, j’ai voulu regarder (inspiré par une remarque de Rémi) si minuit était vraiment une heure de crime1

Ce dont je dispose ce sont les heures d’appels d’urgence (le fameux 911), avec un code (formellement, le code indiqué par la personne lors de l’appel, pas celui établi par la patrouille une fois sur place). Globalement, la distribution des appels (par seconde) au cours de la journée ressemble à ça,

avec les six années en noir, 2005 en bleu et 2010 en rouge (histoire de visualiser une stabilité temporelle).
On peut aussi se focaliser sur les causes données lors de l’appel, afin de voir quel criminel commet un crime à minuit. Par exemple, pour les alarmes consécutives à un cambriolage, on obtient la répartition suivante,

Manifestement, minuit est un creux en matière de cambriolage. N’en déplaise à Arsène Lupin. Les hold-ups, eux aussi surviennent pendant la journée,

Les intrusions-effractions aussi

Sinon, de manière générale, les infractions au code criminel surviennent en fin de journée (mais baissent en fin de soirée)

et les conflits, suivent la même distribution

Autrement dit, les conflits surviennent quand on rentre du travail… mais pas à minuit. En revanche, si on commence à regarder les appels pour signaler descoups de feux, ils surviennent un peu après les conflits… mais un peu avant minuit.

Ah, mais finalement je crois que j’ai trouvé le crime qui survient à minuit,

Et dans ma base, cela correspond au code 063, ce qui signifie trouble à la paix. Autrement dit, si un appel arrive au 911 à minuit, c’est rarement pour signaler un (vrai) crime (au sens de ceux qu’on trouve dans les polars), mais plus pour signaler un problème de beuverie (que l’on peut aussi trouver chez Ian Rankin d’ailleurs). Damned ! Moi qui pensait résoudre des crimes avec mes séries de chiffres…
.
1 Je n’ai pas réussi à trouver d’où pouvait sortir cette idée reçue. En lisant le dernier Ian Rankin (ou plutôt la dernière enquête de John Rebus), je me suis fait la réflexion que dans beaucoup d’enquêtes, les policiers arrivaient vers minuit. J’ai eu l’impression que c’était aussi souvent le cas dans bon nombre d’épisodes de séries télé “policières” (CSI, Castle, Bones, ). En fait si quelqu’un sait où je pourrais trouver une base avec des dates de décès des personnages dans les séries policières ou des polars, je suis preneur.

Will I ever be a bayesian statistician ? (part 2)

A few weeks ago, I started a series of posts on the magic of bayesian statistics from the eyes of a muggle (see http://freakonometrics.hypotheses.org/2191). It might be time to go a bit further…. And today, I wanted to discuss the choice of the a priori distribution of the parameter (which was mentioned in the commentary of the previous post). As far as I understood, there are several houses with different ideas on how to choose it.

  • The conjugate house

The first idea (used in the previous post, here) is to consider an exponential distribution. To be formal, those distributions can be written as

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-10.gif

(in a form as general as possible), i.e.

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-11.gif

Here http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-13.gif is somehow the new parameter of the distribution. Then, a conjugate priorhttp://freakonometrics.blog.free.fr/public/perso2/bayes-prior-12.gif for the parameter http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-13.gif of the exponential family is given by

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-15.gif

where http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-16.gif (where http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-17.gif is the dimension of http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-13.gif).

The conjugate prior is interesting since, when combined with the likelihood (and normalized), produces a posterior distribution which is of the same type as the prior. And a lot of standard distributions have a conjugate prior. E.g.

  • For a Bernoulli distribution, i.e. http://freakonometrics.blog.free.fr/public/perso2/conj-00c.gif are i.i. with distribution http://freakonometrics.blog.free.fr/public/perso2/conju-01.gif, assume that the  a priori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif is http://freakonometrics.blog.free.fr/public/perso2/conju-02.gif, , then the a posteriori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gif is still Beta, with parameters
http://freakonometrics.blog.free.fr/public/perso2/conju-03.gif
  • For a binomial distribution,  http://freakonometrics.blog.free.fr/public/perso2/conju-05.gif, assume that the  a prioridistribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif is http://freakonometrics.blog.free.fr/public/perso2/conju-02.gif, then the a posteriori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gif is still Beta, with parameters
http://freakonometrics.blog.free.fr/public/perso2/conj-06.gif
  • For a Negative Binomial distribution, http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-20.gif,  assume that the  a priori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif is http://freakonometrics.blog.free.fr/public/perso2/conju-02.gif, then the a posterioridistribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gif is still Beta, with parameters
http://freakonometrics.blog.free.fr/public/perso2/conj-08.gif
  • For a Poisson distribution, http://freakonometrics.blog.free.fr/public/perso2/conj-09.gif, assume that the  a priori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif is http://freakonometrics.blog.free.fr/public/perso2/conj-10.gif, then the a posteriori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gifis still gamma, with parameter
http://freakonometrics.blog.free.fr/public/perso2/conj-13.gif
  • For a Geometric distribution,http://freakonometrics.blog.free.fr/public/perso2/conj-15.gif, assume that the  a prioridistribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif is http://freakonometrics.blog.free.fr/public/perso2/conju-02.gif, then the a posteriori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gif is still Beta, with parameters
http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-22.gif
  • For an Exponential distributionhttp://freakonometrics.blog.free.fr/public/perso2/conj-23.gif assume that the  a prioridistribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif ishttp://freakonometrics.blog.free.fr/public/perso2/bayes-prior-21.gif , then the a posteriori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gif is still Gamma, with parameters
http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-22.gif
  • For a Gaussian distribution, http://freakonometrics.blog.free.fr/public/perso2/conj-28.gif assume that http://freakonometrics.blog.free.fr/public/perso2/conj-31.gif, then
http://freakonometrics.blog.free.fr/public/perso2/conj-35.gif

  • For a gamma distribution, http://freakonometrics.blog.free.fr/public/perso2/conj-43.gif, assume that the  a prioridistribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif is  http://freakonometrics.blog.free.fr/public/perso2/bayes-000000.gif, then the a posteriori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gif is still Gamma, with parameters
http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-23.gif
  • For a Pareto distribution, , assume that the  a priori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00b.gif is http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-21.gif, then the a posteriori distribution of http://freakonometrics.blog.free.fr/public/perso2/conj-00.gif is still Gamma, with parameters

http://freakonometrics.blog.free.fr/public/perso2/bayes----ooooo.gif

  • The non-informative or vague house

So far, the choice of the prior was not neutral, in the sense that the a priori of the statistician will have an influence on a posteriori distributions (we’ll discuss that point further later on). We could be interested by the case where http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-12.gif is somehow t neutral. A famous example is the case of  http://freakonometrics.blog.free.fr/public/perso2/conju-01.gifdistribution. Between 1745 and 1784, Pierre Simon Laplace observed 393,386 birth of boys versus 377,555 birth of girls (or 251,527 boys versus 241,945 girls if we consider the initial article, for birth before 1770). He wanted to quantify the probability that p, the provability to have a boy, exceed 1/2. He assume that a priorip was uniform on the unit interval claiming that it was being as neutral as possible. But it is not that correct.
The idea of noninformative prior is that we should get an equivalent result when considering a transformed parameter. So assume that the parameter is no longer http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-03.gif, but http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-05.gif, where http://freakonometrics.blog.free.fr/public/perso2/bayes-----0000.gif (for some bijective transformation). The distribution (density) of http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-05.gif is then

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-01.gif

Let http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-02.gif denote Fisher information of parameter http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-03.gif,

http://freakonometrics.blog.free.fr/public/perso2/bayes-priori-04.gif

Then, Fisher information of parameter http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-05.gif is

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-06.gif

which can be written

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-07.gif

So if we want a distribution invariant by transformation of the parameter (http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-08.gif), it seems natural to consider

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-09.gif

or at least something proportional to that square root, since we want to get a density.
Thus, from Jeffrey’s principle, the prior distribution for a single parameter is noninformative if it is taken proportional to the square root of Fisher’s information measure. For those who want to go further, see Noninformative Priors Do Not Exist or  A Catalog of Noninformative Priors.

  • For the Poisson distribution, the Jeffreys prior for the rate parameter  is
http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-poisson.gif
  • For the Bernoulli distribution, the Jeffreys prior for the probability parameter  is
http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-binomial.gif

This is the arcsine distribution and is a beta distribution with parameters 1/2.

  • The expert house

The idea is quite simple. We need a prior distribution so that

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-32.gif

But assume that we have already seen a similar problem before. For instance, I remember Eric Parent mentioning the case of river flow models. If we have two similar rivers, then it might be interesting to use information on one river as an a prioriinformation for the second one, something like

http://freakonometrics.blog.free.fr/public/perso2/bayes-prior-33.gif

I guess it is also possible to use meta-regression to get an aggregation of experts opinion.

To go further on bayesian statistics, I suggest to go on Albus Dumbledore’s og,here, or the the blog of some PhD (and postdoc) students in Hogwarts, there. Or if you can wait, a dozen other posts will come soon (well, let’s hope so). The next one will probably be on a posteriori calculations (which is the natural step since we’ve seen a priori choice here).

Tests statistiques et intervalles de confiance

Comme je l’avais énoncé rapidement à la fin du chapitre sur les intervalles de confiance, il existe une forme de dualité entre tests et intervalles de confiance. Par exemple, si on considère un test sur la moyenne, de la forme http://freakonometrics.blog.free.fr/public/perso2/testic-01.gifcontre http://freakonometrics.blog.free.fr/public/perso2/testic-02.gif dans un modèle Gaussien, http://freakonometrics.blog.free.fr/public/perso2/test-ic03.gif, la région d’acceptation (de niveau http://freakonometrics.blog.free.fr/public/perso2/testic-04.gif) sera de la forme

http://freakonometrics.blog.free.fr/public/perso2/testic-05.gif
c’est à dire si

http://freakonometrics.blog.free.fr/public/perso2/testic-06.gif
Aussi, pour qu’une valeur http://freakonometrics.blog.free.fr/public/perso2/testic-07.gif soit acceptée avec un niveau http://freakonometrics.blog.free.fr/public/perso2/testic-04.gif (erreur de première espèce associée à un test bilatéral), une condition nécessaire et suffisante est que http://freakonometrics.blog.free.fr/public/perso2/testic-07.gif appartienne à l’intervalle de confiance (symétrique) centré sur la moyenne empirique.
Il y a ainsi dualité entre deux concepts

  • prendre une valeur acceptée pour un test bilatéral de niveau http://freakonometrics.blog.free.fr/public/perso2/testic-04.gif
  • prendre une valeur appartenant à l’intervalle de confiance symétrique de niveau 1-http://freakonometrics.blog.free.fr/public/perso2/testic-04.gif

Considérons par exemple le test asymptotique du rapport de vraisemblance (évoqué ici), dont la région d’acceptation sera de la forme

http://freakonometrics.blog.free.fr/public/perso2/tesp100.gif
Étant donné un échantillon http://freakonometrics.blog.free.fr/public/perso2/testh03.gif, l’intervalle de confiance pour  est

http://freakonometrics.blog.free.fr/public/perso2/tesp101.gif
Numériquement, reprenons l’exemple du jeu de pile ou face abordé ici.

> n=20; alpha=.05
> X=sample(0:1,size=n,replace=TRUE)
> neglogL=function(p){-sum(log(dbinom(X,size=1,prob=p)))}
> pml=optim(fn=neglogL,par=0.5,method="BFGS")$par
Warning messages:
1: In dbinom(x, size, prob, log) : NaNs produced
2: In dbinom(x, size, prob, log) : NaNs produced
> p=seq(0,1,by=.01)
> logL=function(p){sum(log(dbinom(X,size=1,prob=p)))}
>
> TLR=function(p0){2*(logL(pml)-logL(p0))<qchisq(1-alpha,df=1)}
> (IC=sapply(p,TLR))
[1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE
[37]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[49]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[61]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
[73]  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[97] FALSE FALSE FALSE FALSE FALSE
> c(p[min(which(IC==TRUE))],p[max(which(IC==TRUE))])
[1] 0.34 0.75


Graphiquement, on visualise l’intervalle de confiance (à 95%) entre les deux traits verticaux.
De manière duale, à partir d’un intervalle de confiance http://freakonometrics.blog.free.fr/public/perso2/testp-103.gif pour , on peut construire un test. Si http://freakonometrics.blog.free.fr/public/perso2/qqqh20.gif appartient à http://freakonometrics.blog.free.fr/public/perso2/testp-103.gif, on accepte H_0,  sinon on rejette cette hypothèse. Par construction,

http://freakonometrics.blog.free.fr/public/perso2/testp-102.gif
Par exemple l’intervalle de confiance usuel – à 95% – pour une proportion est de la forme suivante

http://freakonometrics.blog.free.fr/public/perso2/testp107.gif
Si http://freakonometrics.blog.free.fr/public/perso2/qqqh20.gif appartient à cet intervalle, on accepte http://freakonometrics.blog.free.fr/public/perso2/testp109.gif (contre l’hypothèse bilatérale http://freakonometrics.blog.free.fr/public/perso2/testp101.gif).

> a=mean(X)-1.96/sqrt(n)*sqrt(mean(X)*(1-mean(X)))
> b=mean(X)+1.96/sqrt(n)*sqrt(mean(X)*(1-mean(X)))
> c(a,b)
[1] 0.3319638 0.7680362


La région entre les deux traits verticaux est la région d’acceptation.

Tests statistiques et p value

Mercredi, nous avons abordé le concept de http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value en cours (avant de parler tout à l’heure de tests asymptotiques). Intuitivement, la http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value est la probabilité que la statistique de test prenne une valeur au moins aussi extrême que celle qui a été observée, si l’hypothèse http://freakonometrics.blog.free.fr/public/perso2/testp-08.gif est vraie.

  • Fisher vs. Neyman (et Pearson)

Dans l’approche de Fisher, on suppose que http://freakonometrics.blog.free.fr/public/perso2/testp-01.gif et on cherche à tester http://freakonometrics.blog.free.fr/public/perso2/TEST-REGION-01.gif. On dispose pour cela d’un échantillon http://freakonometrics.blog.free.fr/public/perso2/testp-02.gif.
On choisit une statistique http://freakonometrics.blog.free.fr/public/perso2/testp-03.gif telle que plus http://freakonometrics.blog.free.fr/public/perso2/testic-5.gif est grande, plus on s’éloigne de http://freakonometrics.blog.free.fr/public/perso2/testp-08.gif (on calcule une distance à l’hypothèse http://freakonometrics.blog.free.fr/public/perso2/testp-08.gif, quelque chose comme ça). On calcule alors la http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value comme

http://freakonometrics.blog.free.fr/public/perso2/testp-07.gif

et on rejette http://freakonometrics.blog.free.fr/public/perso2/testp-08.gif si http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif est trop petite. La notion de “plus petite” est liée à la comparaison avec un seuil prédéfini, souvent 5% (cf plus loin pour une discussion).
C’est un peu heuristique, mais c’est comme cela que ça fonctionne.
Dans l’approche de Neyman-Pearson, on cherche à testerhttp://freakonometrics.blog.free.fr/public/perso2/TEST-REGION-01.gif contre une hypothèse alternative, du genre http://freakonometrics.blog.free.fr/public/perso2/testp-11.gif. Si http://freakonometrics.blog.free.fr/public/perso2/testp-12.gif, on va construire une région de rejet de la forme

http://freakonometrics.blog.free.fr/public/perso2/testp-13.gif

où http://freakonometrics.blog.free.fr/public/perso2/testp-14.gif est choisi en fonction d’une probabilité d’erreur de première espèce souhaitée a priori, i.e.

http://freakonometrics.blog.free.fr/public/perso2/testp-15.gif

Par exemple, considérons un échantillon suivant une loi http://freakonometrics.blog.free.fr/public/perso2/testp-16.gif. On dispose de http://freakonometrics.blog.free.fr/public/perso2/testp-17.gif observations pour tester http://freakonometrics.blog.free.fr/public/perso2/testp-20.gif. L’échantillon est le suivant

> set.seed(4)
> (X=rnorm(10))
 [1]  0.2167549 -0.5424926  0.8911446  0.5959806  1.6356180  0.6892754
 [7] -1.2812466 -0.2131445  1.8965399  1.7768632

La statistique proposée par Fisher (je ne discuterais pas le choix optimal de la statistique – si quelqu’un a des références, les commentaires sont ouvertes) est la moyenne normalisée,

http://freakonometrics.blog.free.fr/public/perso2/testp-21.gif

On sait que sous H_0, cette statistique suit une loi normale centrée réduite.  La http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value se calcule alors simplement

> (T=sqrt(10)*mean(X))
[1] 1.791523
 
> 2*(1-pnorm(T,mean=0,sd=0)))
[1] 0.07320942

Graphiquement, cela se représente de la manière suivante.

Cette probabilité étant faible, on aurait tendance vouloir rejeter http://freakonometrics.blog.free.fr/public/perso2/testp-20.gif.
Dans l’optique de Neyman Pearson, si on teste http://freakonometrics.blog.free.fr/public/perso2/testp-20.gif contre http://freakonometrics.blog.free.fr/public/perso2/testp-22.gif, la région critique s’obtient à l’aide du rapport de vraisemblance

http://freakonometrics.blog.free.fr/public/perso2/testp-25.gif

ou encore

http://freakonometrics.blog.free.fr/public/perso2/testp-26.gif

ce qui se traduit par

http://freakonometrics.blog.free.fr/public/perso2/testp-27.gif

où http://freakonometrics.blog.free.fr/public/perso2/testp-14.gif est tel que

http://freakonometrics.blog.free.fr/public/perso2/testp-29.gif

Or sous H_0, la somme suit une loi normale centrée de variance http://freakonometrics.blog.free.fr/public/perso2/testp-30.gif. Ici, http://freakonometrics.blog.free.fr/public/perso2/testp-14.gif est donné par

> qnorm(1-alpha,mean=0,sd=sqrt(10))
[1] 5.201484

Or pour rappel, la somme vaut ici

> sum(X)
[1] 5.665293

On est dans la région de rejet: on rejette alors http://freakonometrics.blog.free.fr/public/perso2/testp-20.gif.

  • p-value pour les tests usuels de R

La http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value est intéressante en pratique car elle permet à l’utilisateur de ne pas rentrer trop dans les détails de la construction du test, et de la loi de la statistique considérée.
Par exemple (comme celui abordé ici), pour un test d’égalité de proportion, sous R, la commande est

> prop.test(nx,n,.5)
 
1-sample proportions test with continuity correction
 
data:  nx out of n, null probability 0.5
X-squared = 0.05, df = 1, p-value = 0.823
alternative hypothesis: true p is not equal to 0.5
95 percent confidence interval:
0.3204804 0.7617145
sample estimates:
p
0.55

Et de manière générale, tous les tests sous R renvoient une http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value.

  • p ne signifie pas puissance

(mais probabilité). Dans un test simple,http://freakonometrics.blog.free.fr/public/perso2/TEST-REGION-01.gif contrehttp://freakonometrics.blog.free.fr/public/perso2/testh13.gif, la puissance est la probabilité de rejeter http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif sous l’hypothèse où http://freakonometrics.blog.free.fr/public/perso2/H1.gif est vraie (rejeter http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif avec raison si on veut). C’est donc une fonction de  http://freakonometrics.blog.free.fr/public/perso2/testh14.gif.
Pour un test multiple, avec une hypothèse alternative de la forme http://freakonometrics.blog.free.fr/public/perso2/testp-36.gif, on définit alors une fonction puissance pour toute valeur  http://freakonometrics.blog.free.fr/public/perso2/testh14.gif.  On a alors des graphiques de fonction puissance, comme ceux évoqués ici. Bref, la http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value est une valeur, alors que la puissance sera une fonction, définie sur http://freakonometrics.blog.free.fr/public/perso2/testp-40.giflorsque H_1 est une hypothèse de la forme http://freakonometrics.blog.free.fr/public/perso2/testp-41.gif.

  • L’interprétation de p

La http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value est la probabilité d’obtenir (au moins) la statistique obtenue pour l’échantillon si http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif est vraie. Classiquement, on rejette http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif si la http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value est inférieure au seuil http://freakonometrics.blog.free.fr/public/perso2/testh04.gif .
Certains interprètent la http://freakonometrics.blog.free.fr/public/perso2/testic-10.gif-value comme la probabilité que http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif est vraie, http://freakonometrics.blog.free.fr/public/perso2/testp-55.gif. Ce qui n’a pas de sens si l’on s’en tient au formalisme énoncé en cours: http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif est vraie, ou pas.
Ce n’est pas non plus la probabilité de rejeter, à tort, http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif , http://freakonometrics.blog.free.fr/public/perso2/testp-57.gif, correspondant à la puissance du test (cf la discussion dans le paragraphe précédant).

  • Pourquoi 5% ?

En pratique, on prend toujours comme valeur critique 5%, avec la règle décision simple suivante:

  • http://freakonometrics.blog.free.fr/public/perso2/testp-65.gif: si http://freakonometrics.blog.free.fr/public/perso2/testp-61.gif on accepte http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif
  • http://freakonometrics.blog.free.fr/public/perso2/testp-66.gif: si http://freakonometrics.blog.free.fr/public/perso2/testp-60.gif on rejette http://freakonometrics.blog.free.fr/public/perso2/testp-33.gif

C’est Fisher qui semble être à l’origine de ce 5%, en 1925 “The value for which P=0.05, or 1 in 20, is 1.96 or nearly 2; it is convenient to take this point as a limit in judging whether a deviation ought to be considered significant or not. Deviations exceeding twice the standard deviation are thus formally regarded as significant. Using this criterion we should be led to follow up a false indication only once in 22 trials, even if the statistics were the only guide available.” On retrouve ce “un sur vingt” un an plus tard, “Either there is something in the treatment, or a coincidence has occurred such as does not occur more than once in twenty trials”.
Voilà ce que l’on retient souvent de la lecture des conseils de Fisher.
Maintenant, il lui arrive aussi, à l’occasion, de penser que 4% ce n’est pas significatif, e.g. “P is between .02 and .05. The result must be judged significant, though barely so; in view of the data we cannot ignore the possibility that on this field, and in conjunction with the other manures used, nitrate of soda has conserved the fertility better than sulphate of ammonia; the data do not, however, demonstrate this point beyond the possibility of doubt.” Ou parfois pense que 8% est significatif “P=.089. Thus a larger value of 2 would be obtained by chance only 8.9 times in a hundred, from a series of values in random order. There is thus some reason to suspect that the distribution of rainfall in successive years is not wholly fortuitous, but that some slowly changing cause is liable to affect in the same direction the rainfall of a number of consecutive years.”

Wald, score et rapport de vraisemblance

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

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

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

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

ou encore

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

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

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

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

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

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

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

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

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

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

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

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

On pose alors

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

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

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

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

La pente correspond au score,

La statistique de test est alors ici

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

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

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

Playing with quantiles, part 1

A standard idea in extreme value theory (see e.g. here, in French unfortunately) is that to estimate the 99.5% quantile (say), we just need to estimate a quantile of level 95% for observations exceeding the 90% quantile.

In extreme value theory, we assume that the 90% quantile (of the initial distribution) can be obtained easily, e.g. the empirical quantile, and then, for the exceeding observations, we fit a Pareto distribution (a Generalized Pareto one to be precise), and get a parametric quantile for the 95% quantile. I.e.

http://freakonometrics.blog.free.fr/public/perso2/quant01.gif

which can be written

http://freakonometrics.blog.free.fr/public/perso2/quant02.gif

So, an estimation of the cumulative distribution function is

http://freakonometrics.blog.free.fr/public/perso2/quant03.gif

and if we invert it, we get the popular expression for high level quantiles,

http://freakonometrics.blog.free.fr/public/perso2/quant04b.gif

Hence, we do not really care about observations in the core of the distribution.

And I was wondering if this can be transposed with quantile regressions. Hence, I would like to get a quantile regression of level 90% (say) of http://freakonometrics.blog.free.fr/public/perso2/qqq06.gif given http://freakonometrics.blog.free.fr/public/perso2/qqqo5.gif, based on observations http://freakonometrics.blog.free.fr/public/perso2/qqq04.gif‘s, but all observations such that http://freakonometrics.blog.free.fr/public/perso2/qqq07.gif for some http://freakonometrics.blog.free.fr/public/perso2/qqq08.gif are missing. More precisely, I have the following sample (here half of the observations are missing),

Assume that we know that I have observations below the http://freakonometrics.blog.free.fr/public/perso2/qqq06.gif quantile of level 25%, and above the http://freakonometrics.blog.free.fr/public/perso2/qqq06.gif quantile of level 75%.
If I want to get the 90% quantile regression, and the 10% quantile, the code is simply,

library(mnormt)
library(quantreg)
library(splines)
set.seed(1)
mu=c(0,0)
r=0
Sigma <- matrix(c(1,r,r,1), 2, 2)
Z=rmnorm(2500,mu,Sigma)
X=Z[,1]
Y=Z[,2]
 
base=data.frame(X,Y)
plot(X,Y,col="blue",cex=.7)
I=(Y>qnorm(.25)
)&(Y<qnorm(.75))
baseI=base[I==FALSE,]
points(X[I],Y[I],col="light blue",cex=.7)
abline(h=qnorm(.25),lty=2,col="blue")
abline(h=qnorm(.75),lty=2,col="blue")
u=seq(-5,5,by=.02)
reg=rq(Y~X,data=base,tau=.05)
lines(u,predict(reg,newdata=data.frame(X=u)),lty=2)
reg=rq(Y~X,data=baseI,tau=.05*2)
lines(u,predict(reg,newdata=data.frame(X=u)))

The graph is the following

Dotted lines – in black – are theoretical lines (if I had all observations), and plain lines are (where half of the sample if missing). Instead of a standard linear quantile regression, it is also possible to try a spline regression,

So obviously, if I miss something in the middle, that’s no big deal, doted and plain lines are here extremely close.
But what if observations http://freakonometrics.blog.free.fr/public/perso2/qqqo5.gif and http://freakonometrics.blog.free.fr/public/perso2/qqq06.gif were correlated ? Consider a Gaussian random vector http://freakonometrics.blog.free.fr/public/perso2/qqq09.gif with correlation http://freakonometrics.blog.free.fr/public/perso2/qqq10.gif (here 0.
6).

It looks like we overestimate the slope for high quantile, but not for lower quantiles. So if observations are correlated, we have to be cautious with that technique.
But why could that be interesting ? Well, because I wanted to run a quantile regression on marathon results. But I could not get the overall dataset (since I had to import observations manually, and I have to admit that it was a bit boring). So I extracted finish times of the first 10% athletes, and the latest 10%. And I was wondering if it was enough to look at the 5% and 95% quantiles, based on the age of the runner… To be continued.

A Million Random Digits: review of reviews

Recently on his blog (here), Robin mentioned an amazing book, called “A Million Random Digits” published by RAND corporation. The book was initially published in 1955, but RAND published a nice (and expensive) second edition.

A great thing is that on Amazon, there are several extremely interesting reviews of the book. E.g.

4.0 out of 5 stars Didn’t like the ending, February 10, 2009  By Damien Katz

Even though I didn’t really see it coming, the ending was kind of anti-climatic. But overall the book held my attention and I really liked the “10034 56429 234088” part. It’s nice to know I’m not the only one who feels that way.

5.0 out of 5 stars I found a typo, September 14, 2007  By fanfan

To whom do I write to report typographical errors? I noticed that the first “7” on the third line page 48 should be a “3”. The “7” that’s printed there now isn’t random. Other than that, this is really an excellent book.

5.0 out of 5 stars Superb and original plot, April 21, 2007  By Herr Tarquin Biskuitfaß

This one has a very unpredictable plot, sublime character development in a style that stubbornly defies any sort of development in its rare and iconoclastic brilliance, and is told remarkably with numbers instead of letters. Take, for example, this passage on page 202, “98783 24838 39793 80954”. I’m speechless. The symmetry is reminiscent of the I Ching, and it approaches a rare spiritual niveau lacking in American literature. It not only reads well, but it looks great too. I have a tattoo of page 214 on my arm, and I’m hoping to get 202 on my belly to celebrate my next birthday. It is an injustice that Rand Corporation has not received the Nobel Prize for Literature, nor even a Pulitzer.

3.0 out of 5 stars A serious reference work?, October 16, 2006  By BJ

For a supposedly serious reference work the omission of an index is a major impediment. I hope this will be corrected in the next edition.

1.0 out of 5 stars Not Nearly A Million, September 3, 2006  By Liron

This book does not even come close to delivering on its promise of one million random digits. My expectations were high after reading the first sentence, which contained ten unique digits. However, the author seems to have exhasted his creativity in this initial burst, because the other 99.999% of the book is filler in which those same ten digits are shamelessly reused!  If you are looking for a larger offering of numerals in various bases, I highly recommend “Peter Rabbit’s ABC and 123”.

3.0 out of 5 stars Wait for the audiobook version, October 19, 2006  By R. Rosini “Newtype”

While the printed version is good, I would have expected the publisher to have an audiobook version as well. A perfect companion for one’s Ipod.

5.0 out of 5 stars Wait for it…, February 10, 2009  By Cranky Yankee

It started off slow, single digit slow in the beginning but I stuck with it. I eventually learned all about the different numbers, 1,2,3,4,5,6,7,8,9 and 0 and their different combinations.  The author introduced them all a bit too quickly for my taste. I would have been perfectly happy with just 1,2,3,4 and 5 for the first 20,000 digits, but then again, I’m not a famous random-number author, am I?  After a while, patterns emerged and the true nature of the multiverse was revealed to me, and the jokes were kinda funny. I don’t want to spoil anything but you will LOVE the twist ending!  Like 4352204 said to 64231234, “2242 6575 0013 2829!”

Ok, I have to admit I tried to check a few of them (that’s my freaky part). For instance the first one is a fake: the two first numbers – for instance – never show up together (consecutively),

> DIGIT=read.table("
+ http://freakonometrics.blog.free.fr/public/data/digits.txt")
> DIGIT=DIGIT[,2:11]
> k=1
> I=apply(DIGIT[,1:2]==c(10034,56429),1,sum)==2
> for(k in 2:9){
+ I=cbind(I,apply(DIGIT[,k+0:1]==c(10034,56429),1,sum)==2)
+ }
> I0=which(apply(I,1,sum)>0)
> DIGIT[I0,]
 [1] V2  V3  V4  V5  V6  V7  V8  V9  V10 V11
<0 rows> (or 0-length row.names)

Nevertheless, I did have some fun reading those reviews. About the book, unfortunately I have to confess I stopped after 99998 appeared (the first time).