{"word_count": 691, "character_count": 5075}
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,&#62; 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 vraisemblanceou (un peu plus malin) la fonction de log-vraisemblance&#62; p=seq(0,1,by=.01)&#62; logL=function(p)+ {sum(log(dbinom(X,size=1,prob=p)))}&#62; plot(p,Vectorize(logL)(p),+ type="l",col="red",lwd=2)&#62; p0=.5&#62; 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)&#62; neglogL=function(p){-sum(log(dbinom(X,size=1,prob=p)))}&#62; optim(fn=neglogL,par=p0,method="BFGS")$par[1] 0.5499996$value[1] 13.76278Bref 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.&#62; p=NULL&#62; n=20&#62; nx=sum(X==1)&#62; f = expression(nx*log(p)+(n-nx)*log(1-p))&#62; Df = D(f, "p")Si on calcule le score en , comme on s'y attend,&#62; p=.55&#62; eval(Df)[1] 0alors qu'en , on a&#62; p=p0&#62; eval(Df)[1] 4On peut d'ailleurs visualiser l'évolution de cette fonction,&#62; p=seq(0,1,by=.01)&#62; score=eval(Df)&#62; plot(p,score,type="l",col="red",lwd=2)&#62; 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&#62; Df2 = D(Df, "p")&#62; p=p0&#62; -eval(Df2)[1] 80Revenons maintenant aux propriétés du maximum de vraisemblance. Moyennant quelques hypothèses techniques, on a le comportement asymptotique suivantou 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 avoirsoit (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&#62; lancer=matrix(sample(0:1,size=20*1000,+ replace=TRUE),nrow=1000)&#62; xbar=apply(lancer,1,mean)&#62; 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,&#62; lancer=matrix(sample(0:1,size=20*10*10000,+ replace=TRUE),nrow=10000)&#62; xbar=apply(lancer,1,mean)&#62; 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, alorsLe second est lié à un résultat évoqué dans un précédant billet, qui nous permet d'écrire (en supposant que  soit inversible) queEn fait, compte tenu des propriétés de convergence de , on peut aussi écrirePour 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), alorsC'est juste du développement de TaylorA 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 queAvec 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 suivantecar on sait que sous ,Lançons_nous. Le calcul est très simple&#62; Df2 = D(Df, "p")&#62; p=p0&#62; IF=-eval(Df2)&#62; Tobs=(nx/n-p0)*IF*(nx/n-p0)&#62; Tobs[1] 0.2On 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&#62; 1-qchisq(Tobs,df=1)[1] 0.9358152que l'on peut visualiser sur le dessinBon, demain, on poursuit avec le test du rapport de vraisemblance, où on va utiliser la continuité de la log-vraisemblance...