Statistical Tests: Asymptotic, Exact, ou based on Simulations?

This morning, in our mathematical statistics course, we’ve been discussing the ‘proportion test‘, i.e. given a sample of Bernoulli trials$\boldsymbol{x}=\{x_1,\cdots,x_n\}$, with $X_i\sim\mathcal{B}(p)$, we want to test

$H_0:p=p_0$against $H_1:p\neq p_0$

A natural test (which can be related to the maximum likelihood ratio test) is  based on the statistic

$T(\boldsymbol{x})=\sqrt n\frac{\widehat{p} - p_0}{\sqrt{p_0 (1-p_0)}}$

The test function is here

$\psi(\boldsymbol{x})=\boldsymbol{1}(T(\boldsymbol{x})\notin [c_{1,\alpha},c_{2,\alpha}])$To get the bounds of the acceptance region, we need the distribution of $T(\boldsymbol{X})$, under $H_0$. Consider here a numerical application

n=20
p=.5
set.seed(1)
echantillon=sample(0:1,size=n,
prob=c(1-p,p),
replace=TRUE)
• the asymptotic distribution

The first (and standard idea) is to use the central limit theorem, since

$\sqrt n\frac{\hat{p} - p}{\sqrt{p (1-p)}}\overset{\mathcal{L}}{\rightarrow}\mathcal{N}(0,1)$

So, under $H_0$,

$\sqrt n\frac{\hat{p} - p_0}{\sqrt{p_0 (1-p_0)}}\overset{\mathcal{L}}{\rightarrow}\mathcal{N}(0,1)$

Then $c_{1,\alpha}=\Phi^{-1}(\alpha/2)$ while $c_{2,\alpha}=\Phi^{-1}(1-\alpha/2)$. The acceptance region is then between the two red lines, below,

T=sqrt(n)*(mean(echantillon)-.5)/
sqrt(mean(echantillon)*
(1-mean(echantillon)))
u=seq(-3,3,by=.01)
v=dnorm(u)
plot(u,v,type="l",lwd=2)
abline(v=qnorm(.025),col="red")
abline(v=qnorm(.975),col="red")
abline(v=T,col="blue")

• the exact distribution

Here we use the fact that

$\sum_{i=1}^n X_i \sim \mathcal{B}(n,p)$

Using transformation of the ‘density’, we can (at least numerically) compute the (exact) distribution of

$T(\boldsymbol{x})=\sqrt n\frac{\widehat{p} - p_0}{\sqrt{p_0 (1-p_0)}}$

u=seq(-3,3,by=.01)
v=sqrt(.5*(1-.5))*n*dbinom(round(
(sqrt(.5*(1-.5))*u/sqrt(n)+.5)*n),
size=n,prob=.5)/sqrt(n)

Here I used a round value, it guess it would be better with a floor function, but here the graph looks symmetric (which is something I like)

abline(v=sqrt(n)*(qbinom(.025,size=n,prob=.5)/n-.5)/sqrt(.5*(1-.5)),col="red")
abline(v=sqrt(n)*(qbinom(.975,size=n,prob=.5)/n-.5)/sqrt(.5*(1-.5)),col="red")
lines(u,v,type="s")

• distribution based on Monte Carlo simulations

Probably more interesting, here we do not use the fact that we might know the distribution of the mean. We just generate random samples, under $H_0$, and then compute $T(\boldsymbol{X})$,

T=rep(NA,1000)
for(i in 1:1000){
x=sample(0:1,size=n,
prob=c(1-.5,.5),
replace=TRUE)
m=mean(x)
T[i]=(m-.5)/sqrt(m*(1-m))*sqrt(n)}
lines(density(T),lwd=2)
abline(v=quantile(T,.025),col="red")
abline(v=quantile(T,.975),col="red")

Wald, score et rapport de vraisemblance

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

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

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

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

ou encore

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

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

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

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

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

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

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

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

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

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

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

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

On pose alors

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

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

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

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

La pente correspond au score,

La statistique de test est alors ici

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

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

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