Freakonometrics

To content | To menu | To search

Teaching › statistique 10/11 STT2700

Entries feed - Comments feed

Thursday, July 14 2011

Comparer des notes entre pays

A Montréal, l'an passé, le principal souvenir que je garderai de mon cours est le mal que j'ai eu à comprendre comment mettre des notes. Je ne ferais pas un billet sur la difficulté à évaluer le niveau d'un étudiant (car j'ai toujours détesté les examens, aussi bien en tant qu’étudiant que prof), mais j'ai une une réelle difficulté à comprendre ce que vaut une note... Par exemple, un C, ça correspond à quoi (avec mon référentiel de prof français) ? C'est bien ou pas ? C'est quoi la moyenne ?
La seule information dont on dispose pour mettre les notes est le tableau ci-contre. Et autant je pense voir ce qu'est un travail "excellent", ou "faible", mais "passable" ? Car le tableau me dit qu'un élève qui a 50 sur 100 a une note passable, ce qui est aussi la référence que l'on utilise en France. Sauf que les élèves ont commence à me dire que 50, c'est pas terrible... qu'une note moyenne ici c’était plutôt 60, voire 65. Bref, sous-entendu (voire entendu, tout court), "monsieur, vous devez remonter les notes". Et comme je l'ai prouve ici, il ne faut jamais croire les élèves...
Et effectivement, dans un vieux rapport datant des années 90 (en maths et informatiques), je suis tombé sur le graphique suivant

qui donne la distribution "typique" des notes. Autrement dit, la médiane est à 63, ce qui signifie que plus de moitie a une note qui est "bien","très bien" voire "excellent". En revanche, une distribution comme celle dessous est considérée comme "dure et discriminatoire"

Pour mieux comprendre les notations au Canada et en France, j'ai voulu faire une petite étude (non, je ne reviendrai pas sur les notes). On continue aujourd'hui à trouver ce genre de distribution à Montréal. Par exemple, si on considère les cours, analyse (en première année), algèbre linéaire et calcul, on a les histogrammes suivants

ce qui donne le boxplot suivant

i.e. pour les cours algèbre linéaire et calcul on est effectivement sur des moyennes autour de 60.
Prenons maintenant, a titre de comparaison, des cours de maths, dans une université française, en l’occurrence a Rennes, comme analyse (en première année), arithmétique et algèbre linéaire, on a des histogrammes beaucoup plus plats,

A titre de comparaison, j'ai été voir du cote du département de Sciences Économiques, avec macroéconomie, microéconomie et mathématiques appliquées, on a

i.e. on retrouve de belles gaussiennes pour les cours d’économie, mais une distribution assez différente pour le cours de maths,

On a effectivement des distributions assez différentes, généralement, en première année à l’université, entre le Canada et la France. Ce qui explique que les profs, lorsque les étudiants souhaitent passer d'un pays a l'autre ne regarde plus vraiment les notes, mais plus le rang de l’étudiant au sein de sa promotion.
D'ailleurs, si on fait un graphique quantile-quantile, avec les cours de maths à Montréal versus les cours de maths a Rennes, on a le graphique suivant

autrement dit, une note de 50 à Montréal correspond à 33 en France (soit 6.5 sur 20), alors que 50 en France (10 sur 20) correspond a 59 à Montréal. On retrouve la encore le fait qu'une note médiane correspond à 60 ici... c'est a dire un C. Il va falloir que je m'adapte l’année prochaine....

Wednesday, May 4 2011

STT2700, examen final - éléments de correction

L'examen final du cours de statistique avait lieu en fin de semaine dernière. Le sujet est en ligne ici, et j'ai posté des éléments de correction (il doit probablement resté quelques coquilles ici ou là malgré tout). Les notes seront envoyées dans la journée à la scolarité.

Concernant le dernier devoir maison, comme il a été corrigé en classe le dernier jour, je ne pensais pas mettre en ligne de correction.

Wednesday, April 13 2011

STT2700, examen final

L'examen final pour le cours Concepts et méthodes en statistique (STT2700) aura lieu Vendredi 29 avril, à 13:30 et durera 3 heures (détails ici) salle Z310. L'examen portera sur l'ensemble du cours, avec 1 exercice de probabilité, 3 exercices d'inférence, et 2 exercices sur les tests. Les exercices sont tirés du livre de John Rice, Mathematical Statistics and Data Analysis, Duxbury Press, 3rd ed. [QA 276.12]. Comme pour l'intra, je donnerais les formes des lois, sauf pour les lois usuelles qui sont supposées connues (loi binomiale, normale, et Poisson). Sont supposés connus

  • les calculs d'espérance et de variance de n'importe quelle loi
  • les théorèmes limites (lois des grands nombres et théorème central limite) et les méthodes qui tournent autour (delta method)
  • la construction d'estimateurs et l'obtention de propriétés (biais et comportement asymptotique)
  • les concepts autour des tests (erreurs de première et seconde espèce, puissance), le théorème de Neyman-Pearson et la construction de région de rejet
Les méthodes bayésiennes, les tests asymptotiques (Wald ou score) et les tests nonparamétriques (Kolmogorov Smirnov) ne sont pas au programme.


Monday, April 11 2011

STT2700, Comparer des taux sur deux populations (indépendantes)

Lors du dernier cours, nous avions évoqué les tests de comparaison des moyennes et des proportions sur deux échantillons. Pour cela, nous avions vu qu'il était possible d'utiliser la statistique

http://freakonometrics.blog.free.fr/public/perso2/comp-sample.gif

http://freakonometrics.blog.free.fr/public/perso2/compar-sample2.gif

La statistique doit alors suivre, sous l'hypothèse où les proportions sont égales une loi normale centrée réduite. Sous R, c'est assez facile à implémenter. Afin d'illustrer, utiliser l'information suivante

Bon, en l'état on n'a pas vraiment un taux, mais un nombre de tués sur les routes. On peut introduire la probabilité d'avoir un tué sur les routes par minute. Le code sous R serait alors
> prop.test(x=c(308,308/1.027), n=c(31*24*60,31*24*60), alt="two.sided")
 
2-sample test for equality of proportions with continuity correction
 
data: c(308, 308/1.027) out of c(31 * 24 * 60, 31 * 24 * 60)
X-squared = 0.0834, df = 1, p-value = 0.7727
alternative hypothesis: two.sided
95 percent confidence interval:
-0.0009198487 0.0012826342
sample estimates:
prop 1 prop 2
0.006899642 0.006718249
Autrement dit, la probabilité d'avoir un accident, à la minute, est sensiblement le même, avec une p-value de 77%. On notera sur le graphique suivant que l'impact du pas de temps a globalement peu d'importance....)

L'idéal serait d'avoir des données détaillées, mais il faudra attendre un peu pour ça...

Friday, April 1 2011

STT2700, estimation, tests et coupes du monde

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 http://freakonometrics.free.fr/test-soccer-01.gif. On va pouvoir tester
http://freakonometrics.free.fr/test-soccer-04.gif
contre une hypothèse alternative
http://freakonometrics.free.fr/test-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 http://freakonometrics.free.fr/test-soccer-03.gif, 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
http://freakonometrics.free.fr/test-soccer-11.gif
Un test de l'hypothèse

http://freakonometrics.free.fr/test-soccer-04.gif
contre une hypothèse alternative
http://freakonometrics.free.fr/test-soccer-06.gif
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

Tuesday, March 29 2011

STT2700, ajustement de loi, test de Kolmogorov-Smirnov, PP et QQ plots

Pour finir le chapitre sur les tests statistiques, on avait abordé la problématique de l'ajustement de loi, qui est un test non plus sur un paramètre (dans http://freakonometrics.free.fr/ks-gif-01.gif voire http://freakonometrics.free.fr/ks-gif-02.gif), mais sur les fonctions de répartition (dans le cas simple http://freakonometrics.free.fr/ks-gif-03.gif). Formellement, on quitte le monde de la statistique paramétrique pour entrer dans celui de la statistique non-paramétrique. Afin de poursuivre ce que j'ai fait en cours, je voulais montrer comment obtenir, numériquement, la loi limite qui apparait dans le test de Kolomogorov Smirnov.
  • La fonction de répartition empirique
La fonction de répartition obtenue est partir d'un échantillon i.i.d. est
http://freakonometrics.blog.free.fr/public/maths/GlivC01.gif
Sous R, la fonction est
> ecdf(X)
Empirical CDF
Call: ecdf(X)
x[1:25] = -1.2846, -1.2375, -0.89192, ..., 1.0858, 1.1519
  • Un mot sur les PP ou QQ plots
Un outil (graphique) particulièrement intéressant pour faire des tests d'ajustement de loi sont les PP (probabilité versus probabilité) plot, et les QQ (quantile versus quantile) plot. L'idée est simple. On se demande si, pour tout http://freakonometrics.free.fr/ks-gif-04.gif, F(x) est proche de http://freakonometrics.free.fr/ks-gif-05.gif. Étant donné que http://freakonometrics.free.fr/ks-gif-06.gif "converge" (on reviendra sur une formalisation de ce concept dans 2 minutes) vers http://freakonometrics.free.fr/ks-gif-07.gif, on va se demander si, pour tout http://freakonometrics.free.fr/ks-gif-04.gif, http://freakonometrics.free.fr/ks-gif-08.gif est proche de http://freakonometrics.free.fr/ks-gif-05.gif. Au lieu de se poser la question pour tout http://freakonometrics.free.fr/ks-gif-09.gif, on peut simplement se demander si pour tout http://freakonometrics.free.fr/ks-gif-10.gif http://freakonometrics.free.fr/ks-gif-11.gif est proche de http://freakonometrics.free.fr/ks-gif-12.gif. Ou de manière équivalente, si la loi http://freakonometrics.free.fr/ks-gif-13.gif est continue, si http://freakonometrics.free.fr/ks-gif-14.gif est proche de http://freakonometrics.free.fr/ks-gif-15.gif.
On peut alors faire un graphique représentant le nuage de points
http://freakonometrics.free.fr/ks-gif-16.gif
c'est à dire, tout simplement
http://freakonometrics.free.fr/ks-gif-17.gif
compte tenu de ce que nous avions noté dans le paragraphe précédant. Comme on représente des probabilités (les fonctions de répartition sont des probabilités), on parlera de PP-plot.
X=rnorm(n)
Xs=sort(X)
plot(pnorm(Xs,0,1),(1:n)/n)
Si les points sont alignés suivant la première diagonale, on acceptera l'hypothèse que http://freakonometrics.free.fr/ks-gif-18.gif est la loi des observations.
On peut toutefois lire le test non pas sur les fonctions de répartition mais sur les quantiles. On va alors se demander si, pour tout http://freakonometrics.free.fr/ks-gif-19.gif, http://freakonometrics.free.fr/ks-gif-20.gif est proche de http://freakonometrics.free.fr/ks-gif-34.gif. Étant donné que http://freakonometrics.free.fr/ks-gif-06.gif "converge" vers http://freakonometrics.free.fr/ks-gif-07.gif, et donc que http://freakonometrics.free.fr/ks-gif-30.gif "converge" vers http://freakonometrics.free.fr/ks-gif-31.gif, on aura que, pour tout http://freakonometrics.free.fr/ks-gif-19.gif, http://freakonometrics.free.fr/ks-gif-33.gif sera proche de http://freakonometrics.free.fr/ks-gif-34.gif. Au lieu de se poser la question pour tout u, on peut simplement se demander si pour tout http://freakonometrics.free.fr/ks-gif-10.gif, http://freakonometrics.free.fr/ks-gif-35.gif est proche de http://freakonometrics.free.fr/ks-gif-37.gif. Ou de manière équivalente, si http://freakonometrics.free.fr/ks-gif-38.gif est proche de http://freakonometrics.free.fr/ks-gif-39.gif.
On peut alors faire un graphique représentant le nuage de points
http://freakonometrics.free.fr/ks-gif-40.gif
Comme on représente des quantiles, on parlera de QQ-plot.
plot(qnorm((1:n)/n,0,1),Xs)
Là encore, si les points sont alignés suivant la première diagonale, on acceptera l'hypothèse que http://freakonometrics.free.fr/ks-gif-18.gif est la loi des observations.

Il existe toutefois quelques cas particulier. Comme le cas Gaussien. En effet, si http://freakonometrics.free.fr/ks-gif-18.gif est la loi http://freakonometrics.free.fr/ks-gif-41.gif, alors
http://freakonometrics.free.fr/ks-gif-42.gif
http://freakonometrics.free.fr/ks-gif-52.gif est la fonction de répartition de la loi normale centrée réduite. On peut alors, dans le cas Gaussien (ou supposé tel), représenter le nuage de points
http://freakonometrics.free.fr/ks-gif-55.gif
Si les points sont alignés suivant une droite, de pente http://freakonometrics.free.fr/ks-gif-57.gif, passant par http://freakonometrics.free.fr/ks-gif-56.gif à l'origine, alors http://freakonometrics.free.fr/ks-gif-18.gif sera la loi http://freakonometrics.free.fr/ks-gif-41.gif.
Et plus généralement, si les points sont alignés, alors on acceptera la normalité des observations.
Remarque: cette propriété sur la loi normale se retrouve aussi sur la loi exponentielle, par exemple.
Mais ce test est graphique... On peut avoir envie d'aller un peu plus loin, et d'introduire un test plus formel permettant de savoir quand la distance entre le nuage de points et la droite est trop grande, ou au contraire suffisamment petite pour valider l'ajustement.

  • Théorèmes limites, à http://freakonometrics.free.fr/ks-gif-61.gif donné, i.e. sur http://freakonometrics.free.fr/ks-gif-60.gif
A http://freakonometrics.free.fr/ks-gif-09.gif donné, on a un modèle binomial, au sens où http://freakonometrics.blog.free.fr/public/maths/binom_cv_01.gif suit une loi http://freakonometrics.blog.free.fr/public/maths/binom_cv_02.gif, et donc
http://freakonometrics.blog.free.fr/public/maths/binom_cv_03.gif
La loi des grands nombres garantie que http://freakonometrics.blog.free.fr/public/maths/binom_cv_05.gif et le théorème central limite que
http://freakonometrics.blog.free.fr/public/maths/binom_cv_04.gif
  • Théorèmes limites, sur distances entre fonctions, i.e. sur http://freakonometrics.free.fr/ks-gif-65.gif
Si la loi est continue, on peut montrer un analogue de la loi de grands nombres, au sens où
http://freakonometrics.blog.free.fr/public/maths/GlivC02.gif
En normalisant, on peut obtenir une convergence en loi, i.e. une sorte de théorème centrale limite,
http://freakonometrics.blog.free.fr/public/maths/GlivC09.gif
On notera que l'on peut translater le problème, par un changement de variable bijectif (en utilisant la fonction quantile)
http://freakonometrics.blog.free.fr/public/maths/GlivC08.gif
  • La loi limite dans le test de Kolmogorov-Smirnov
Mais il est plus simple de faire des simulations pour trouver la loi limite (et trouver les seuils de régions critiques des tests),
http://freakonometrics.free.fr/animationks.gif

> n=25
> D=rep(NA,10000)
> Utheo1=(0:(n-1))/n
> Utheo2=(1:n)/n
> for(s in 1:10000){
+ Uemp=sort(runif(n))
+ D[s]=max(c(abs(Uemp-Utheo1),abs(Uemp-Utheo2)))
+ }
> quantile(D,.95)
95%
0.2644086
On peut ainsi récupérer numériquement une bonne approximation d'une p-value (la valeur exacte est obtenue à l'aide du test de R),
> set.seed(1)
> ks.test(rnorm(n),pnorm)
 
One-sample Kolmogorov-Smirnov test
 
data: rnorm(n)
D = 0.2021, p-value = 0.2263
alternative hypothesis: two-sided
 
 
> mean(D>0.2021)
[1] 0.2309

STT2700, changement de variables

Afin d'illustrer une des questions du devoir maison (en ligne ici), je vais revenir sur le changement de variables multivarié (évoqué très rapidement ). Je l'avais fait en cours dans le cas univarié, et je pensais que rappeler la formule dans le cas multivarié suffirait. Mais on va faire un exemple pratique (pas celui demandé dans le devoir maison, mais l'idée est la même), ça ne fera pas de mal...

  • Retour sur la théorie

Soit http://freakonometrics.free.fr/chg-gif-01.gif un couple de variables aléatoires continues (ou plutôt absolument continues), de telle sorte que la loi du couple admet une densité de probabilité http://freakonometrics.free.fr/chg-gif-02.gif .  Soit http://freakonometrics.free.fr/chg-gif-03.gif un changement de variables bijectif. Si on note http://freakonometrics.free.fr/chg-gif-04.gif le Jacobien associé à la transformation, i.e.

http://freakonometrics.free.fr/chg-gif-05.gif
alors la loi de http://freakonometrics.free.fr/chg-gif-06.gif est donnée par la densité
http://freakonometrics.free.fr/chg-gif-09.gif
Ca c'est pour la théorie.
  • Application(s) aux statistiques d'ordre 
Soient http://freakonometrics.free.fr/chg-gif-10.gif une collection de http://freakonometrics.free.fr/chg-gif-12.gif variables aléatoires indépendantes, positives, de même loi de densité http://freakonometrics.free.fr/chg-gif-11.gif.Soit http://freakonometrics.free.fr/chg-gif-13.gif le réarrangement par ordre croissant des variables http://freakonometrics.free.fr/chg-gif-10.gif, i.e.
http://freakonometrics.free.fr/chg-gif-14.gif
Notons http://freakonometrics.free.fr/chg-gif-15.gif et notons http://freakonometrics.free.fr/chg-gif-16.gif la différence entre deux statistiques d'ordre avec la convention http://freakonometrics.free.fr/chg-gif-17.gif. Le vecteur http://freakonometrics.free.fr/chg-gif-18.gif peut alors être vu comme un réarrangement du vecteur http://freakonometrics.free.fr/chg-gif-19.gif pour un certaine permutation http://freakonometrics.free.fr/chg-gif-20.gif de http://freakonometrics.free.fr/chg-gif-21.gif. Le nombre de permutation de http://freakonometrics.free.fr/chg-gif-21.gif étant http://freakonometrics.free.fr/chg-gif-22.gif, on en déduit que
http://freakonometrics.free.fr/chg-gif-23.gif
Pour détailler un peu, cette relation peut se retrouver en considérant http://freakonometrics.free.fr/chg-gif-12.gif=2 par exemple. Considérons  http://freakonometrics.free.fr/chg-gif-24.gif, et alors, en notant abusivement http://freakonometrics.free.fr/chg-gif-25.gif (on va se placer dans le cas discret pour illustrer). Alors
http://freakonometrics.free.fr/chg-----02.gif
soit tout simplement
http://freakonometrics.free.fr/chg-gif-27.gif
Bref, la loi de la statistique d'ordre est simple à obtenir, et finalement assez intuitive. Mais maintenant on veut la loi des distances entre statistiques d'ordre. Pour faire simple, plaçons nous dans le cas où les http://freakonometrics.free.fr/chg-gif-28.gif sont distribués suivant une loi exponentielle, alors http://freakonometrics.free.fr/chg-gif-30.gif, et donc, la loi de la statistique d'ordre est
http://freakonometrics.free.fr/chg-gif-31.gif
A partir de cette densité, il est possible d'obtenir la loi du vecteur http://freakonometrics.free.fr/chg-gif-32.gif. On peut pour cela utiliser la formule de changement de variable évoquée dans le paragraphe précédant. Ici,
http://freakonometrics.free.fr/chg-gif-40.gif
qui peut s'inverser en
http://freakonometrics.free.fr/chg-gif-41.gif
Le Jacobien vaut alors 1,
http://freakonometrics.free.fr/chg-gif-45.gif
Pour l'ensemble dans l'indicatrice, http://freakonometrics.free.fr/chg-gif-42.gif, notons qu'il se réécrit http://freakonometrics.free.fr/chg-gif-43.gif et la densité devient
http://freakonometrics.free.fr/chg----05.gif

http://freakonometrics.free.fr/chg----04.gif
On en déduit (par la propriété de séparabilité) que les variables http://freakonometrics.free.fr/chg-gif-55.gif sont indépendantes, et que leur loi est une loi exponentielle. On en déduit alors par exemple que
http://freakonometrics.free.fr/chg-gif-60.gif
et que
http://freakonometrics.free.fr/chg-gif-63.gif

Friday, March 25 2011

STT2700, 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.

Tuesday, March 22 2011

STT2700, 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.blog.free.fr/public/maths/chi2-16.gif, et on cherchait à tester une hypothèse de la forme
http://freakonometrics.blog.free.fr/public/maths/chi2-13.gif
versus
http://freakonometrics.blog.free.fr/public/maths/chi2-14.gif
Ici on va cherche à utiliser un test sur une loi multinomiale, de la forme
http://freakonometrics.blog.free.fr/public/maths/chi2-01.gif
versus
http://freakonometrics.blog.free.fr/public/maths/chi2-02.gif
L'hypothèse nulle est ici une égalité vectorielle,
http://freakonometrics.blog.free.fr/public/maths/chi2-03.gif
ou encore
http://freakonometrics.blog.free.fr/public/maths/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.blog.free.fr/public/maths/chi2-09.gif, on va utiliser un test sur une loi multinomiale, avec
http://freakonometrics.blog.free.fr/public/maths/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.blog.free.fr/public/maths/chi2-10.gif
Le test est alors basé sur la statistique de Pearson
http://freakonometrics.blog.free.fr/public/maths/chi2-11.gif
qui va suivre (si effectivement les observations suivent une loi de Poisson http://freakonometrics.blog.free.fr/public/maths/chi2-09.gif) une loi http://freakonometrics.blog.free.fr/public/maths/chi2-12.gif.
Rappelons que cette statistique peut aussi s'écrire
http://freakonometrics.blog.free.fr/public/maths/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.blog.free.fr/public/perso2/Nb.gif le nombre de points tombés dans un ensemble http://freakonometrics.blog.free.fr/public/perso2/asubb.gif (ou http://freakonometrics.blog.free.fr/public/perso2/calA.gif désigne la région globale). Si on suppose qu'un nombre aléatoire http://freakonometrics.blog.free.fr/public/perso2/Npois.gif de points sont lancés aléatoirement dans http://freakonometrics.blog.free.fr/public/perso2/calA.gif, et que http://freakonometrics.blog.free.fr/public/perso2/Npois.gif suit une loi de Poisson http://freakonometrics.blog.free.fr/public/perso2/lambda.gif, alors http://freakonometrics.blog.free.fr/public/perso2/Nb.gif suit une loi de Poisson de paramètre
http://freakonometrics.blog.free.fr/public/perso2/ratio-aire.gif
Bref, si on regarde plusieurs plusieurs régions http://freakonometrics.blog.free.fr/public/perso2/calB.gif (de même taille, éventuellement - afin de garder toutes les observations - formant une partition de http://freakonometrics.blog.free.fr/public/perso2/calA.gif) et si on observe une loi de Poisson, c'est que dans la région http://freakonometrics.blog.free.fr/public/perso2/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.blog.free.fr/public/perso2/.logvrais1-bombes_m.jpg
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.blog.free.fr/public/perso2/.chi2-bombes_m.jpg
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( ((freq.theo-freq.emp)^2)/freq.theo ))
[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( ((n/sum(n)-c(dpois(y,L),1-ppois(4,L)))^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.blog.free.fr/public/perso2/.chi2-bombes-v2_m.jpg
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...

Saturday, March 19 2011

STT2700, 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://freakonometrics.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.

Thursday, March 17 2011

Statistiques, STT2700, Devoir Maison No 2

Le deuxième devoir maison sera à rendre pour vendredi 1er avril 2011, l'énoncé est en ligne ici (le jour du dernier TD). Quelques billets seront postés afin de donner des lignes de codes permettant de calculer - numériquement - certaines quantités.

Pour importer la base, le code est le suivant

x=c(.6,.19,6.44,1.74,.02,2.34,4.08,.17,
1.16,1.23,1.74,.55,.25,3.43,1.64,5.32,
3.8,1.27,.68,2.22,2.77,2.85,3.67,.26,
2.48,4.08,1.23,.35,5.05,3.22,0)
sum(x)
[1] 64.83
sum(x^2)
[1] 227.1353

(merci Maxim qui a noté que la somme des carrés était fausse.)

Thursday, March 10 2011

Statistiques, STT2700, 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.gif contre 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.

Statistiques, STT2700, 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
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
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.gif lorsque 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."

Statistiques, STT2700, 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 . Alors
(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.

Sunday, March 6 2011

Statistiques, STT2700, introduction aux tests

En cours, avant la semaine de relâche, on a commencé la théorie des tests. Dans le cours, nous avons introduit les tests comme une prise de décision. On se demande si une hypothèse - notée http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif - est vraie, ou pas. Et plus précisément, on note http://freakonometrics.blog.free.fr/public/perso2/H1.gif l'hypothèse dite alternative, correspondant à l'hypothèse qui doit être vérifiée sur http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif n'est pas vérifiée.
A partir d'un échantillon http://freakonometrics.blog.free.fr/public/perso2/testh03.gif, on doit prendre une décision. Schématiquement, on peut prendre la décision http://freakonometrics.blog.free.fr/public/perso2/test-D0.gif, i.e. accepter l'hypothèse http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif, ou bien la décision http://freakonometrics.blog.free.fr/public/perso2/test-D1.gif, i.e. on accepter l'hypothèse http://freakonometrics.blog.free.fr/public/perso2/H1.gif (ou rejeter http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif).
La décision (d'accepter ou de rejeter une hypothèse) est faite en fonction d'une région où peut se trouver - ou pas - l'échantillon. Schématiquement, un test consiste à construire une région d'acceptation (ou plutôt de rejet), i.e.

http://freakonometrics.blog.free.fr/public/perso2/test-region-1.gif
http://freakonometrics.blog.free.fr/public/perso2/test-intro-2.gif
De manière équivalente, on peut chercher à construire une statistique telle que
http://freakonometrics.blog.free.fr/public/perso2/test-intro-4.gif
http://freakonometrics.blog.free.fr/public/perso2/test-intro-5.gif
Dans les tests paramétriques, un exemple classique de tests est
http://freakonometrics.blog.free.fr/public/perso2/TEST-REGION-01.gif
avec comme hypothèse alternative
http://freakonometrics.blog.free.fr/public/perso2/TEST-REGION-2.gif
(mais on le verra, plusieurs autres formes sont possibles).
Une fois construit un test (i.e. un mécanisme de prise de décision), il faut se demander si le test fonctionne bien. Pour cela, on peut noter que quatre situations peuvent exister,

http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif http://freakonometrics.blog.free.fr/public/perso2/H1.gif
http://freakonometrics.blog.free.fr/public/perso2/test-D0.gif bonne décision
erreur de seconde espèce
http://freakonometrics.blog.free.fr/public/perso2/test-D1.gif erreur de première espèce
bonne décision
Dans deux cas, la décision est bonne, mais dans deux cas, on a commis une erreur.
Le premier type d'erreur consiste à rejeter notre hypothèse http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif alors qu'elle était valide. On note alors http://freakonometrics.blog.free.fr/public/perso2/testh04.gif la probabilité de faire une erreur de première espèce,
http://freakonometrics.blog.free.fr/public/perso2/test-intro-7.gif
Le second type d'erreur consiste à accepter notre hypothèse http://freakonometrics.blog.free.fr/public/perso2/test-H0.gif à tort, et la probabilité faire une erreur de seconde espèce est http://freakonometrics.blog.free.fr/public/perso2/testh05.gif
http://freakonometrics.blog.free.fr/public/perso2/test-intro-8.gif
Lors que l'on construit un test, on va fixer l'erreur de premier type, et regarder ce que peut être l'erreur de second type.
Dans le cas d'un test  sur la moyenne d'un échantillon normal http://freakonometrics.blog.free.fr/public/perso2/testh06.gif, et que l'on cherche à tester
http://freakonometrics.blog.free.fr/public/perso2/TEST-REGION-01.gif
avec comme hypothèse alternative
http://freakonometrics.blog.free.fr/public/perso2/TEST-REGION-2.gif
La statistique (usuelle) de ce test est
http://freakonometrics.blog.free.fr/public/perso2/testh07.gif
Comme on l'avait noté il y a quelques semaines, si H_0 est vraie,
http://freakonometrics.blog.free.fr/public/perso2/testh08.gif
La région de rejet (la meilleur, comme le montre le théorème de Neyman Pearson) de niveau \alpha est de la forme suivante

i.e. la région critique (ou de rejet) est
http://freakonometrics.blog.free.fr/public/perso2/testh09.gif
http://freakonometrics.blog.free.fr/public/perso2/testh10.gif est le quantile de la loi de Student à http://freakonometrics.blog.free.fr/public/perso2/testh11.gif degrés de liberté, de niveau http://freakonometrics.blog.free.fr/public/perso2/testh12.gif.
Pour étudier l'erreur de seconde espèce associée à ce test, on va se considérer une hypothèse alternative simple
http://freakonometrics.blog.free.fr/public/perso2/testh13.gif
La probabilité d'erreur de seconde espèce est l'aire jaune ci-dessous

pour une valeur dehttp://freakonometrics.blog.free.fr/public/perso2/testh14.gif proche de http://freakonometrics.blog.free.fr/public/perso2/qqqh20.gif, ou pour une valeur relativement différente,

En fait, au lieu de calculer l'erreur de seconde espèce, on s'intéresse à la puissance définie par
puissancehttp://freakonometrics.blog.free.fr/public/perso2/testh15.gif
(que l'on souhaite la plus grande possible). Le graphique ci-dessous montre la puissance de ce test avec http://freakonometrics.blog.free.fr/public/perso2/testh16.gif observations

en fonction de la valeur de la différence entre http://freakonometrics.blog.free.fr/public/perso2/testh14.gif et http://freakonometrics.blog.free.fr/public/perso2/qqqh20.gif (en abscisse). Plus les valeurs sont différentes, meilleure sera le test. Ce n'est pas pour nous déplaire.
Voilà pour la terminologie générale de la théorie de la décision vue par Neyman et Pearson. Je reviendrais bientôt sur Fisher et la p value. à suivre...

- page 1 of 3