Petit rappel préliminaire. Si on dispose d’un échantillon , i.id de loi , où le paramètre  est inconnu, alors on peut calculer le maximum de la vraisemblance, ce qui nous donnera un estimateur intéressant (cf. cours de stat). Illustrons avec un jeu de pile ou face. Reprenons l’échantillon de mon précédent billet,

> X=c(0, 0, 1, 1, 0, 1, 1, 1, 1, 0,
      0, 0, 1, 0, 1, 0, 1, 1, 0, 1)

On peut ici tracer la fonction de vraisemblance

ou (un peu plus malin) la fonction de log-vraisemblance

> p=seq(0,1,by=.01)
> logL=function(p)
+ {sum(log(dbinom(X,size=1,prob=p)))}
> plot(p,Vectorize(logL)(p),
+ type="l",col="red",lwd=2)
> p0=.5
> abline(v=p0,col="blue")

La valeur correspondant au trait bleu correspondant au cas que l’on va chercher à tester, de pièce non pipée, soit  (de manière très générale)

Si on va chercher le maximum de la vraisemblance, on obtient (il faut changer de signe car R a tendance à chercher des minimums, et pas des maximums)

> neglogL=function(p){-sum(log(dbinom(X,size=1,prob=p)))}
> optim(fn=neglogL,par=p0,method="BFGS")
$par
[1] 0.5499996

$value
[1] 13.76278

Bref le maximum de la vraisemblance semble être atteint quand  vaut 0.55. Comme par hasard, on notera que c’est la moyenne de mon échantillon, et effectivement, dans le cas d’une variable de Bernoulli, le maximum de vraisemblance est l’estimateur des moments.

Si on va un peu plus loin, on a deux quantités importantes quand on manipule la fonction de vraisemblance. La première est le score, qui est la pente de la log-vraisemblance,

On peut utiliser un propriété de R, qui peut – à l’occasion – faire un peu de calcul formel.

> p=NULL
> n=20
> nx=sum(X==1)
> f = expression(nx*log(p)+(n-nx)*log(1-p))
> Df = D(f, "p")

Si on calcule le score en , comme on s’y attend,

> p=.55
> eval(Df)
[1] 0

alors qu’en , on a

> p=p0
> eval(Df)
[1] 4

On peut d’ailleurs visualiser l’évolution de cette fonction,

> p=seq(0,1,by=.01)
> score=eval(Df)
> plot(p,score,type="l",col="red",lwd=2)
> abline(v=p0,col="blue")

L’autre quantité presque plus importante, c’est la dérivée seconde (ou la Hessienne en dimension plus grande)

Par exemple en ,

prend la valeur

> Df2 = D(Df, "p")
> p=p0
> -eval(Df2)
[1] 80

Revenons maintenant aux propriétés du maximum de vraisemblance. Moyennant quelques hypothèses techniques, on a le comportement asymptotique suivant

ou encore, en posant  (cf. cours sur l’information de Fisher, c’est d’ailleurs cette dernière valeur que l’on a calculée, numériquement)

Sous l’hypothèse , on devrait donc avoir

soit (puisqu’on a ici un paramètre univarié)

D’un point de vue heuristique, cela signifie que si je pouvais lancer plusieurs fois des pièces que je sais non-pipées, je devrais pouvoir générer un échantillon de loi normale centrée et réduite. Essayons pour voir

> lancer=matrix(sample(0:1,size=20*1000,
+ replace=TRUE),nrow=1000)
> xbar=apply(lancer,1,mean)
> hist(sqrt(800)*(xbar-p0),
+ breaks=seq(-4,4,by=.5),probability=TRUE)

Bon, je l’admets, ce n’est pas super Gaussien… mais en fait, c’est juste que la normalité est asymptotique, et ici, on fait seulement 20 lancers. Si on en fait 10 fois plus,

> lancer=matrix(sample(0:1,size=20*10*10000,
+ replace=TRUE),nrow=10000)
> xbar=apply(lancer,1,mean)
> hist(sqrt(80*10)*(xbar-p0),
+ breaks=seq(-4,4,by=.5),probability=TRUE)

c’est mieux, non ?

Maintenant, allons un peu plus loin pour présenter les trois tests classiques construits autour de la vraisemblance. Mais deux résultats importants avant d’aller plus loin. Le premier est la delta-methode, qui repose juste sur un développement de Taylor. L’idée est que si on a une fonction  suffisemment régulière, alors

Le second est lié à un résultat évoqué dans un précédant billet, qui nous permet d’écrire (en supposant que  soit inversible) que

En fait, compte tenu des propriétés de convergence de , on peut aussi écrire

Pour des paramètres unidimensionnels, on dit juste que si

(avec une notation plus standard car dans les notations précédantes, la taille de l’échantillon était cachée dans la matrice d’information de Fisher), alors

C’est juste du développement de Taylor

A condition toutefois que . Sinon, il faut aller un cran plus loin dans le développement,

et dans ce cas, au lieu d’avoir une loi normale, à droite, on a une loi du chi-deux,

de telle sorte que

Avec ces propriétés, on peut enfin faire des tests. En particulier pour tester des hypothèses simple, du genre

 contre  

Le premier test que l’on  a vu est le test de Wald.

L’idée est d’utiliser la statistique de test suivante

car on sait que sous ,

Lançons_nous. Le calcul est très simple

> Df2 = D(Df, "p")
> p=p0
> IF=-eval(Df2)
> Tobs=(nx/n-p0)*IF*(nx/n-p0)
> Tobs
[1] 0.2

On peut alors calculer la p-value, en calculant la probabilité qu’une loi du chi-deux à un degré de libeté (on a un paramètre unidimensionnel) dépasse cette value

> 1-qchisq(Tobs,df=1)
[1] 0.9358152

que l’on peut visualiser sur le dessin

Bon, demain, on poursuit avec le test du rapport de vraisemblance, où on va utiliser la continuité de la log-vraisemblance…

---

---