Tag Archives: chi-deux

Test du chi-deux et indépendance

22/10/2014 Arthur Charpentier Leave a comment

Considérons le tableau de contingence suivant

> N=margin.table(HairEyeColor, c(1,2))
> N
       Eye
Hair    Brown Blue Hazel Green
  Black    68   20    15     5
  Brown   119   84    54    29
  Red      26   17    14    14
  Blond     7   94    10    16

avec ici les comptages, que l’on peut aussi traduire sous la forme de probabilités jointes (empiriques)

> N/n
       Eye
Hair    Brown  Blue Hazel Green
  Black 0.115 0.034 0.025 0.008
  Brown 0.201 0.142 0.091 0.049
  Red   0.044 0.029 0.024 0.024
  Blond 0.012 0.159 0.017 0.027

Continue reading Test du chi-deux et indépendance →

M1-statistique, Statistics, Tests

Loi multinomiale et loi du chi-deux

15/10/2014 Arthur Charpentier Leave a comment

La semaine passée, en cours, j’avais rappelé que quand $https://latex.codecogs.com/gif.latex?\boldsymbol{N}=(N_{1},\cdots,N_{k})$ décrivait le compte de $https://latex.codecogs.com/gif.latex?n$ variable multinomiales prenant $https://latex.codecogs.com/gif.latex?k$ modalités, la variable

$https://latex.codecogs.com/gif.latex?\tilde{X}_i%20=%20\frac{N_i-np_i}{\sqrt{np_i\color{red}{(1-p_i)}}}$

suit asymptotiquement une loi $https://latex.codecogs.com/gif.latex?\mathcal{N}(0,1)$ . Et plus généralement, on peut montrer que $https://latex.codecogs.com/gif.latex?\tilde{\boldsymbol{X}}\sim\mathcal{N}(\boldsymbol{0},\star)$ . Le soucis est que la matrice de variance covariance n’est pas la matrice identité. Pire que ça, elle n’est pas diagonale. Encore pire, elle n’est pas inversible.

On ne peut alors pas utiliser le joli résultat qui nous dit que si $https://latex.codecogs.com/gif.latex?\boldsymbol{X}\sim\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ , alors $https://latex.codecogs.com/gif.latex?Y=(\boldsymbol{X}-\boldsymbol{\mu})^{\text{\sffamily%20T}}\boldsymbol{\Sigma}^{-1}(\boldsymbol{X}-\boldsymbol{\mu})$ suit une loi du chi-deux à $https://latex.codecogs.com/gif.latex?k$ degrés de liberté (dimension du vecteur $https://latex.codecogs.com/gif.latex?\boldsymbol{X}$ ). Car il faut que $https://latex.codecogs.com/gif.latex?\boldsymbol{\Sigma}$ soit inversible. Or ici, je pense qu’on peut montrer que

$https://latex.codecogs.com/gif.latex?\det(\boldsymbol{\Sigma})=\prod_{i=1}^k%20p_i%20\cdot\left(\sum_{i=1}^k%20p_i%20-1\right)$

(à un facteur $https://latex.codecogs.com/gif.latex?n$ près). Et comme on a une loi de probabilité, le terme de droite est nul. Bref, $https://latex.codecogs.com/gif.latex?\det(\boldsymbol{\Sigma})=0$ et son inverse n’existe alors pas. Cela dit, $https://latex.codecogs.com/gif.latex?0$ n’est valeur propre qu’une fois, d’où la dimension $https://latex.codecogs.com/gif.latex?k-1$ qui apparaît sans cesse (à condition que les $https://latex.codecogs.com/gif.latex?p_i$ soient non-nuls, ce qui était une hypothèse dans la définition de la loi multinomiale : toutes les modalités existent vraiment).

L’idée est alors de poser

$https://latex.codecogs.com/gif.latex?X_i%20=%20\frac{N_i-np_i}{\sqrt{np_i}}$

suit asymptotiquement une loi $https://latex.codecogs.com/gif.latex?\mathcal{N}(0,\star)$ . Et plus généralement, là encore, $https://latex.codecogs.com/gif.latex?\boldsymbol{X}\sim\mathcal{N}(\boldsymbol{0},\star)$ .

Ce qui est amusant, c’est que dans le cas binomiale (j’avais dit qu’il était intéressant de voir la loi multinomiale comme une généralisation de la loi binomiale),

$https://latex.codecogs.com/gif.latex?\tilde{X}_1^2=X_1^2+X_2^2$

puisque

$https://latex.codecogs.com/gif.latex?\frac{(N_1-np_1)^2}{np_1(1-p_1)}=\frac{(N_1-np_1)^2}{np_1}+\frac{(N_2-np_2)^2}{np_2}$

Et en fait, un peut de calculs rapides permettent de noter que

$https://latex.codecogs.com/gif.latex?(N_1-np_1)^2=(N_2-np_2)^2$

de telle sorte que $https://latex.codecogs.com/gif.latex?\tilde{X}_1^2=\tilde{X}_2^2$ alors que $https://latex.codecogs.com/gif.latex?{X}_1^2\neq%20{X}_2^2$ . Bref, c’est clairement sur cette dernière quantité qu’il faut travailler !

Et là on fait des calculs.

L’astuce est d’utiliser la formule de Stirling

$https://latex.codecogs.com/gif.latex?n!\sim%20n^{n+\frac{1}{2}}e^{-n}\sqrt{2\pi}$

alors pour $https://latex.codecogs.com/gif.latex?n_1+\cdots+n_k=n$

$https://latex.codecogs.com/gif.latex?\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})=\frac{n!}{n_1%20!%20\cdots%20n_k%20!}p_1^{n_1}\cdots%20p_k^{n_k}$ soit

$https://latex.codecogs.com/gif.latex?\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})\sim\frac{n^{n+\frac{1}{2}}e^{-n}\sqrt{2\pi}}{\prod_{j=1}^k%20[n_j^{n_j+\frac{1}{2}}e^{-n_j}\sqrt{2\pi}]}p_1^{n_1}\cdots%20p_k^{n_k}$

et donc

$https://latex.codecogs.com/gif.latex?\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})\sim%20\kappa%20\prod_{j=1}^k%20\left(\frac{np_j}{n_j}\right)^{n_j+\frac{1}{2}}$

de telle sorte que

$https://latex.codecogs.com/gif.latex?\log%20\mathbb{P}(\boldsymbol{N}=\boldsymbol{n})=k+\sum_{j=1}^k%20\left(n_j+\frac{1}{2}\right)%20\log%20%20\left(\frac{np_j}{n_j}\right)$ On peut alors faire des développements limités, en faisant intervenir

$https://latex.codecogs.com/gif.latex?x_i=\frac{n_i-np_i}{\sqrt{np_i}}$

et non pas la version avec $https://latex.codecogs.com/gif.latex?\sqrt{np_i(1-p_i)}$ au dénominateur. On peut alors réécrire la somme en faisant intervenir $https://latex.codecogs.com/gif.latex?k-1$ termes (et non pas $https://latex.codecogs.com/gif.latex?k$ ), exactement comme dans les calculs faits dans le premier cours, avec les calculs de variances. Bref, un peu de calcul supplémentaire permet de montrer que

$https://latex.codecogs.com/gif.latex?\sum_{i=1}^k%20X_i^2=\sum_{i=1}^k%20\frac{(N_i-np_i)^2}{np_i}\overset{\mathcal{L}}{\rightarrow}%20\chi^2(k-1)$

C’est donc

$https://latex.codecogs.com/gif.latex?Q=\sum_{i=1}^k%20\frac{(N_i-np_i)^2}{np_i}$

et non pas

$https://latex.codecogs.com/gif.latex?\tilde{Q}=\sum_{i=1}^k%20\frac{(N_i-np_i)^2}{np_i\color{red}{(1-p_i)}}$

qui sera utilisée comme statistique de test par la suite.

Pour utiliser le test du chi-deux comme test d’ajustement de loi, je pourrais reprendre l’exemple des buts lors des coupes du monde de Football

soccer1982=read.table("http://freakonometrics.free.fr/soccer1982")
S82=(soccer1982$V1+soccer1982$V2)
soccer1998=read.table("http://freakonometrics.free.fr/soccer1998")
S98=(soccer1998$V1+soccer1998$V2)
soccer2010=read.table("http://freakonometrics.free.fr/soccer2010")
S10=(soccer2010$V1+soccer2010$V3)

ACT6420, Regression

sur les transformations dans un modèle linéaire

16/11/2012 Arthur Charpentier 1 Comment

Je voulais prendre 5 minutes pour reprendre une question posée par courriel, qui me permettra de poursuivre sur des choses évoquées mercredi dernier en cours. Je vais reformuler la question, mais en gros, cela disait: la méthode de Box-Cox (évoquée en début de semaine, ici) avait pour objet de choisir entre un modèle linéaire (que l’on étudié dans tous les sens depuis le début du cours) et un modèle log-linéaire (évoqué ici, par exemple). L’idée était d’associer le cas linéaire à $\lambda=1$ (dans la transformée de Box-Cox) et à $\lambda=0$ pour le cas multiplicatif (ou log-linéaire). Mais que se passe-t-il si la valeur optimale rejette ces deux cas, et est proche (disons) de $1/2$ ?
Considérons le cas suivant (exemple que je ne me lasse d’utiliser)

> reglm=lm(dist~speed,data=cars)
> library(MASS)
> boxcox(reglm)

La valeur optimale est effectivement proche de $1/2$ . En posant $Z_i=Y_i^{1/2}$ , on aurait envie de considérer un modèle de la forme

$Z_i=\sqrt{Y_i}=\beta_0+\beta_1 X_i+\varepsilon_i$

(comme le suggère la transformation de Box-Cox). On peut faire la régression,

> regsqrt=lm(sqrt(dist)~speed,data=cars)
+ summary(regsqrt)

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

Residuals:
Min      1Q  Median      3Q     Max
-2.0684 -0.6983 -0.1799  0.5909  3.1534

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept)  1.27705    0.48444   2.636   0.0113 *
speed        0.32241    0.02978  10.825 1.77e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.102 on 48 degrees of freedom
Multiple R-squared: 0.7094,	Adjusted R-squared: 0.7034
F-statistic: 117.2 on 1 and 48 DF,  p-value: 1.773e-14

et effectivement, le résultat est concluant. Mais on ne peut pas en rester là, comme pour pour le passage au logarithme, notre variable d’intéret reste $Y$ .
La prédiction sur le modèle transformé était

$\widehat{Z}_i=\widehat{\beta}_0+\widehat{\beta_1}X_i$

$\mathbb{E}(\widehat{Z}_i)=Z_i$

Comme avec le logarithme, on pourrait prendre comme prédiction pour $\widehat{Y}_i=\widehat{Z}_i^2$ mais

$\mathbb{E}(\widehat{Y}_i)=\mathbb{E}(\widehat{Z}_i^2)\neq \mathbb{E}(\widehat{Z}_i)^2=Z_i^2=Y_i$

Bref, on a un estimateur biaisé. Et là encore, comme on est sur des variables positives, l’inégalité de Jensen doit meme nous garantir que l’une des quantités domine toujours l’autre (je laisse les amateurs de convexité l’écrire dans le bon sens). Comment faire ?
La première solution est de noter que

$\mathbb{E}(Z^2)=\text{Var}(Z)+\mathbb{E}(Z)^2$

Donc de manière très simple, on a notre prédiction (en prenant le carré, comme intuité, mais en rajoutant une terme – positif – lié à la variance).
Une seconde consiste à noter que Z (conditionnellement à la variable explicative) était gaussienne. Donc en prenant le carré (moyennant quelques changement d’échelle), on devrait tomber sur une loi du chi-deux qui est une loi que l’on connait bien (sinon je peux renvoyer ici).
Essayons de creuser un peu ces idées (en particulier, pour dériver des intervalles de confiance). Si on visualise la prédiction de ce modèle sur la racine carrée de la distance de freinage, on obtient

> plot(speed,sqrt(dist))
> x=seq(0,30,by=.2)
> distsqrtp=predict(regsqrt,newdata=
+ data.frame(speed=x),interval="prediction")
> polygon(c(x,rev(x)),c(distsqrtp[,2],
+ rev(distsqrtp[,3])),col="yellow",border=NA)
> lines(x,distsqrtp[,1],lwd=2,col="red")

Pour passer d’une loi normale à une loi du chi-deux, il faut retrancher la moyenne (pour centrer la variable) et diviser par l’écart-type (pour avoir une variance unitaire),

> s=summary(regsqrt)$sigma
> mu=predict(regsqrt)
> distsqrtp01=(sqrt(dist)-mu)/s

On a ainsi une variable qui, conditionnellement à la variable explicative est supposé suivre une loi $\mathcal{N}(0,1)$ ,

> plot(speed,distsqrtp01,ylim=c(-3,3))
> abline(h=qnorm(c(.025,.975),0,1),lty=2,col="red")

On prend ensuite le carré pour obtenir notre loi $\chi^2(1)$ , on on trace la bande de confiance (il faut que la valeur à la deux soit faible),

> plot(speed,distsqrtp01^2)
> abline(h=qchisq(.95,df=1),lty=2,col="red")

On a maitenant nos intervalles de confiance, en utilisant cette loi du chi-deux, en inversant notre transformation,

> mu=predict(regsqrt,newdata=data.frame(speed=x))
> distsup=(mu+s*sqrt(qchisq(.95,df=1)))^2
> distinf=(mu-s*sqrt(qchisq(.95,df=1)))^2

que l’on peut visualiser simplement (sur les données de base, avec en ordonnée la distance)

> plot(cars)
> lines(x,distsup,lty=2,col="red")
> lines(x,distinf,lty=2,col="red")

Et la valeur prédite ? On va utiliser notre relation liant la variance, l’espérance du carré, et le carré de l’espérance,

> distesp=mu^2+s^2
> lines(x,distesp,lwd=2,col="red")

Nice, n’est-ce pas ?
Sinon, comme je le disais en cours, si la transformation optimale sur $Y$ est de prendre $\lambda=1/2$ , peut-etre pourrait-on envisager de prendre $X^2$ (de manière un peu duale) comme variable explicative. Encore une fois, c’est une idée, car rien ne le garantit. En effet

$(\beta_0+\beta_1X+\varepsilon)^2 \neq \alpha_0+\alpha_1X+\alpha_2X^2+\eta$

En particulier, le second modèle est un modèle linéaire classique, avec de l’homoscédasticité: on aura une dispersion uniforme autour de notre parabole. En revanche, avec le premier modèle, on a une espèce d’hétéroscédasticité (comme tenu du double produit), avec une variance du terme d’erreur qui croit avec $X$ (car les coefficients sont positifs). Bref, en terme d’intervalles de confiance, on devrait avoir des choses assez différentes.
Regardons la régression de la distance (cette fois) sur la vitesse, et le carré de la vitesse,

> reglm=lm(dist~speed+I(speed^2),data=cars)
> distp=predict(reglm,newdata=
+ data.frame(speed=x),interval="prediction")

Si on visualise la prédiction, avec un intervalle de confiance, on obtient

> plot(cars)
> polygon(c(x,rev(x)),c(distp[,2],
+ rev(distp[,3])),col="yellow",border=NA)
> lines(x,distp[,1],lwd=2,col="red")

Autrement dit, on suppose que notre modèle est homoscédastique. Les régions de confiance entre les deux approches sont clairement différentes… meme si les prédictions sont presque superposées (oui oui, il y a deux courbes, une rouge et une bleue).

Computer, Statistics, STT2700

STT2700, estimation, tests et coupes du monde

01/04/2011 Arthur Charpentier Leave a comment

Mercredi, pour le dernier cours, nous allons revenir sur l’estimation, les tests, et plus généralement sur la modélisation statistique. Pour cela, j’avais pensé travailler sur les nombres de buts marqués, par match, lors de différentes coupes du monde de soccer (1982, 1998 et 2010). Je ne mets pas l’intégralité du code aujourd’hui, l’idée est pour l’instant de mettre en ligne des données qui serviront à répondre aux questions qui seront posées mercredi. Le code (accompagné – éventuellement – d’explications théoriques) sera posté par la suite.

soccer1982=read.table("http://freakonometrics.free.fr/soccer1982")
S82=(soccer1982$V1+soccer1982$V2)
soccer1998=read.table("http://freakonometrics.free.fr/soccer1998")
S98=(soccer1998$V1+soccer1998$V2)
soccer2010=read.table("http://freakonometrics.free.fr/soccer2010")
S10=(soccer2010$V1+soccer2010$V3)

Les boxplot associés à ces trois échantillons sont les suivants,

On va se poser des questions autour de ces données, par exemple voir s’il est vraisemblance (ou pas) que le nombre moyen de but dans un match (avant prolongation, s’il y en a eu). On peut commencé par essayer de se demander quel modèle utiliser. Classiquement, la loi de Poisson est la plus utilisée (en plus, c’est la seule loi qui est autorisée lorsqu’on publie un billet le 1er avril). Les histogrammes sont les suivants

boxplot(S82,S98,S10,col=c("red","yellow","blue"),
label=("1982","1998","2010"))
hist(S82,breaks=0:11,col="red")
hist(S98,breaks=0:11,col="yellow")
hist(S10,breaks=0:11,col="blue")

Si on compare les fonctions de répartition empiriques à celles de lois de Poisson ajustées par maximum de vraisemblance, on obtient, pour la coupe du monde de 1982

et pour celle de 2010,

Visuellement, l’ajustement semble relativement bon, surtout en 2010. On peut aussi faire un test du chi-deux,

> library(vcd)
> (GF=goodfit(S10,type="poisson"))

Observed and fitted values for poisson distribution
with parameters estimated by ML

 count observed     fitted
     0        7  6.6409703
     1       17 15.0459484
     2       13 17.0442384
     3       14 12.8719508
     4        7  7.2907534
     5        5  3.3036226
     6        0  1.2474617
     7        1  0.4037543

> summary(GF)

	 Goodness-of-fit test for poisson distribution

                      X^2 df  P(> X^2)
Likelihood Ratio 5.586765  5 0.3485255

On voit que l’on accepte l’ajustement par une loi de Poisson. Pour ceux qui veulent une visualisation, sur la figure ci-dessous, on a la densité d’une loi du chi-deux. Le premier trait vertical est la valeur observée, et l’aire jaune est alors la p-value (qui excède largement 5%). En rouge on a 5%, donc le second trait vertical est la borne de la région critique associé au test pour une erreur de première espèce valant 5%,

On peut ensuite faire plein de tests. On suppose que $X_i\sim\mathcal{P}(\lambda)$ . On va pouvoir tester

http://freakonometrics.free.fr/test-soccer-04.gif

contre une hypothèse alternative

http://freakonometrics.free.fr/test-<br /><br /> soccer-06.gif

Comme on a une hypothèse sur la loi des observations qui semble robuste, on peut utiliser un test de type rapport de vraisemblance.
On peut aussi faire un test de la forme

http://freakonometrics.free.fr/test-soccer-09.gif

contre

http://freakonometrics.free.fr/test-soccer-10.gif

(histoire de tester des hypothèses simples – qui ont une interprétation). Sinon, comme ce qui nous intéresse, c’est de savoir si on a plus de trois buts par matchs, on peut définir la variable binomiale , en notant que

http://freakonometrics.free.fr/test-soccer-02.gif

est une proportion – donc facile à tester – qui nous intéresse ici compte tenu du problème que l’on cherchera à résoudre. On pourra alors tester, par exemple

http://freakonometrics.free.fr/test-soccer-08.gif

contre

http://freakonometrics.free.fr/test-soccer-07.gif

Ces derniers tests sont alors facile à mettre en œuvre,

> Z=(S10>=3)*1
> prop.test(sum(Z),length(Z),p=1/2,alternative="less")

	1-sample proportions test with continuity correction

data:  sum(Z) out of length(Z), null probability 1/2 
X-squared = 1.2656, df = 1, p-value = 0.1303
alternative hypothesis: true p is less than 0.5 
95 percent confidence interval:
 0.0000000 0.5322764 
sample estimates:
       p 
0.421875

On peut aussi faire des tests de moyenne sur

Un test de l’hypothèse

contre une hypothèse alternative

s’écrit alors

> t.test(S10,mu=3,alternative ="less")

	One Sample t-test

data:  S10 
t = -3.7763, df = 63, p-value = 0.0001775
alternative hypothesis: true mean is less than 3 
95 percent confidence interval:
     -Inf 2.590273 
sample estimates:
mean of x 
 2.265625

Mais on l’aura l’occasion de revoir tous les points du cours, y compris aller peut être un peu plus loin, par exemple sur la comparaison de moyenne entre échantillons,

> t.test(S82,S98,var.equal=FALSE)

	Welch Two Sample t-test

data:  S82 and S98 
t = 0.427, df = 85.266, p-value = 0.6704
alternative hypothesis: true difference in means is not equal to 0 
95 percent confidence interval:
 -0.5503669  0.8514658 
sample estimates:
mean of x mean of y 
 2.807692  2.657143

Inference, Statistics, STT2700

Test d’ajustement et lâcher de bombes

22/03/2011 Arthur Charpentier 1 Comment

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 , 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 , 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 ) une loi .
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 km² 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 km² 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 le nombre de points tombés dans un ensemble (ou désigne la région globale). Si on suppose qu’un nombre aléatoire de points sont lancés aléatoirement dans , et que suit une loi de Poisson , alors 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 (de même taille, éventuellement – afin de garder toutes les observations – formant une partition de ) et si on observe une loi de Poisson, c’est que dans la région , 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.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)/
+ 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…

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

Computer, Inference, STT2700, Tests

Cochrane, Pearson et le test du chi-deux

19/03/2011 Arthur Charpentier 2 Comments

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 peut prendre modalités, notées , avec . On posera

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

en notant que appartient au simplexe de 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 ). Mais plus généralement, on finira les calculs permettant d’établir que, pour

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

(ce qui permet d’obtenir la matrice de variance covariance de ).
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

où

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 est .
Aussi, 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 est une matrice de projection si elle est idempotente, i.e. . Ses valeurs propres sont alors 0 ou 1, et si est le nombre de fois où 1 est valeur propre, et si , alors .
Posons . 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-40.gif

qui est une matrice de projection dont la trace est (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 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 de dimension . Posons , pour , où on notera le rang de , en supposant que les sont positives semidefinies, alors on a équivalence entre

pour ,
les sont des variables indépendantes.

Les s’interprètent comme des longueurs (euclidienne) de projections d’un vecteur Gaussien sur des sous-espaces orthogonaux (de dimension respective ). Si ces vecteurs sont indépendants, et suivent des lois du chi-deux à degrés de libertés, avec , 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 pouvant prendre toutes les deux deux modalités (disons deux lois binomiales). On est alors face a une loi multinomiale à 4 modalités

avec probabilité
avec probabilité
avec probabilité
avec probabilité

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 et . tels que et . On a alors 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 i.e. 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.

Climate

Des étés pluvieux en Bretagne ? une réalité statistique…

21/08/2010 Arthur Charpentier Leave a comment

Pour compléter le précédant billet (ici) on peut se demander en quoi la Bretagne est différente des autres régions françaises… Nous avions vu ici le niveau de précipitation moyen, jour après jours pendant les mois d’été, en Bretagne. A Rennes. En revanche, à Paris on obtient la moyenne suivante,

que l’on peut comparer à Marseille,

ou encore à Strasbourg,

Sur la figure ci-dessous, on voit que la probabilité d’avoir de la pluie à Paris (au sens au moins 0.1 mm d’eau dans la journée, en trait gras bleu, au moins 2 mm d’eau dans la journée, en trait bleu) est supérieure à la probabilité d’avoir de la pluie à Rennes (respectivement en bleu clair gras, et en bleu clair fin)

On est certes très au dessus de Marseille,

mais très en dessous de Strasbourg,

Mais au delà des lois marginales, ces villes sont différentes de la Bretagne si l’on regarde les matrices de transition.

Transition d’un jour sur l’autre

Pour Rennes, si on regarde jour après jour, on obtient

		jour
		beau temps	pluie
jour	beau temps	1955	612	2567
jour	pluie	606	723	1329