{"word_count": 706, "character_count": 5848}
Jeudi, dans le cadre du cours de statistique mathématiques, nous avions passé un peu de temps à mettre en oeuvre le test de Wald, et discuter la construction et l'interprétation d'un test. On reste ici sur un test paramètrique. La vraie valeur du paramètre, on va la fixer, pour générer des données&#62; p=.5En l'occurence, on va considérer un modèle de Bernoulli, i.e. une succession de tirages de pièces (pile ou face)&#62; n=20&#62; set.seed(1)&#62; echantillon=sample(c("pile","face"),+    size=n,replace=TRUE,+    prob=c(p,1-p))&#62; echantillon[1] "face" "face" "pile" "pile" "face" "pile" "pile" "pile" "pile" "face" "face" "face" "pile" "face" "pile" "face" "pile" "pile" "face" "pile"Classiquement, on va recoder en &#62; X=(echantillon=="pile")*1&#62; X[1] 0 0 1 1 0 1 1 1 1 0 0 0 1 0 1 0 1 1 0 1La vraisemblance est ici donnée par&#62; L=function(p) prod(dbinom(X,1,p))&#62; logL=function(p) sum(log(dbinom(X,1,p)))&#62; xp=seq(0,1,by=.01)&#62; yL=Vectorize(L)(xp)&#62; plot(xp,yL,type="l")mais on va plutôt travailler sur la log-vraisemblance&#62; ylogL=Vectorize(logL)(xp)&#62; plot(xp,ylogL,type="l")Le maximum de cette fonction est obtenu, numériquement, avec&#62; optim(par=.2,fn=function(p) -logL(p))$par[1] 0.55$value[1] 13.76278Le maximum est ici atteint en 0,55, ce qui correspond d'ailleurs à l'estiimateur de la méthode des moments&#62; mean(X)[1] 0.55Passons maintenant à la phase de test. Notre hypothèse  est que la pièce n'est pas truquée,&#62; p0 = 1/2     #  Hypoth?se H0&#62; phat = mean(X)Commençons par calculer l'information de Fisher&#62; IF=function(p) (p*(1-p))^(-1)&#62; IF(phat)[1] 4.040404On a utilisé l'expression vue en cours, car on a fait le calcul, au tableau. Mais on peut aussi faire un peu de calcul formel (et utiliser le faire que la dérivée seconde est linéaire, de telle sorte que l'espérance de la dérivée seconde de la log vraisemblance est la valeur de la dérivée seconde en l'espérance)&#62; nx=sum(X==1)&#62; f = expression(nx*log(p)+(n-nx)*log(1-p))&#62; df = D(f,"p")&#62; d2f = D(df,"p")&#62; p = phat&#62; -eval(d2f)[1] 80.80808ce qui correspond bien à l'information de Fisher pour l'échantillon&#62; n*IF(phat)[1] 80.80808On peut alors calculer la statistique de Wald. La procédure de test est basée sur le quantile d'une loi du chi-deux à un degré de liberté (qui correspond à la loi - asymtptotique - sous  de la statistique de Wald)&#62; T=n*IF(phat)*(phat-p0)^2&#62; seuil=qchisq(p=.95,df=1)&#62; T&#60;=seuil # TRUE : accepter H0[1] TRUEOn accepte l'hypothèse   si la statistique est inférieure au seuil. Une autre solution est de calculer la p value&#62; pvalue=1-pchisq(T,df=1)&#62; pvalue&#62;0.05 # TRUE : accepter H0[1] TRUEOn rejette l'hypothèse   si la p-value est inférieure à 5%. 5% est ici le seuil de significativité, correspondant à la probabilité de commettre une erreur de première espèce, i.e. rejeter à tort l'hypothèse  . Pour interpréter ce 5%, on peut générer des échantillons, sous   (i.e. avec la probabilité 1/2), et compter la proportion des scénarios où la statistique de test dépasse le seuil&#62; simu_Wald=function(graine){+   set.seed(graine)+   X=sample(0:1,size=n,replace=TRUE)+   phat=mean(X)+   T=n*IF(phat)*(phat-1/2)^2+   decision=T&#62;qchisq(p=.95,df=1)+   pvalue=1-pchisq(T,df=1)+   return(list(+     echantillon=X,+     param=phat,+     stat=T,+     decision=decision*1,+     pvalue=pvalue))+ }&#62;&#62; D=Vectorize(function(s){+   simu_Wald(s)$decision})(1:1e5)&#62; mean(D)[1] 0.04094On rejette à tort   dans environ 5% des cas. Ce qui correspond effectivement au choix du seuil que nous avions fait.On peut aussi visualier les valeurs des moyennes obtenues dans les 1000 premiers échantillons, et les représenter en rouge quand on a pris la décision de rejeter &#62; P=Vectorize(function(s){+   simu_Wald(s)$param})(1:1e5)&#62; plot(P[1:1000],1:1000,col=D[1:1000]+1)On s'aperçoit ici que le test de Wald est un test sur la moyenne. Ici, si la proportion de 'pile' est comprise entre 0,3 et 0,7, on va accepter l'hypothèse  .On peut aussi regarder la distribution de la p-value&#62; PV=Vectorize(function(s){+   simu_Wald(s)$pvalue})(1:1e5)&#62; hist(PV)On est loin de la loi uniforme ici, mais ça vient du fait qu'on est sur une loi discrète (avec peu d'observations). Cela dit, la p-value excède 5% dans 95% des cas...&#62; mean(PV&#62;.05)[1] 0.95906On peut ensuite regarder la distribution de la statistique de test&#62; St=Vectorize(function(s){+   simu_Wald(s)$stat})(1:1e5)&#62; hist(St,prob=TRUE, breaks =seq(0,50,by=.2),+      xlim=c(0,10))&#62; x=seq(0,50,by=.2)&#62; lines(x,dchisq(x,df=1),lwd=2,col="red")Là encore, on est loin de la loi du chi-deux, mais c'est normal: c'est asymptotiquement que la statistique de Wald suit une loi du chi-deux. Notons juste que le quantile empirique, à 95% est ici&#62; quantile(St,.95)95%3.809524qui n'est pas loi du quantile théorique de la loi du chi-deux&#62; seuil[1] 3.841459Finallement, on peut regarder la (fonction) puissance de ce test&#62; simu_puissance=function(graine,p1){+   set.seed(graine)+   X=sample(0:1,size=n,replace=TRUE,+            prob=c(1-p1,p1))+   phat=mean(X)+   T=n*IF(phat)*(phat-1/2)^2+   decision=T&#62;qchisq(p=.95,df=1)+   return(decision)+ }&#62; puissance=function(p1){+   p=Vectorize(function(s){+     simu_puissance(s,p1)})(1:1e5)+   mean(p)+ }&#62; P1=seq(.05,.95,by=.05)&#62; PUISSANCE=Vectorize(puissance)(P1)&#62; plot(P1,PUISSANCE,type="b")où on retrouve l'allure classique de la puissance d'un test bilatéral.Pour conclure, notons que la statistique de Wald est (numériquement) proche de celle du test du rapport de vraisemblance&#62; T=2*(logL(phat)-logL(p0))&#62; T[1] 0.2003347et du test du score&#62; sc=function(x,p) x/p-(1-x)/(1-p)&#62; Sc=Vectorize(function(x) sc(x,p0))(X)&#62; T=n*(mean(Sc))^2/IF(phat)&#62; T[1] 0.198ce qui est rassurant...