Un cadre comparatif basé sur une « expérience de pensée »

Afin de comparer des méthodes d'estimations, les statisticiens ont recours à ce qu'on désignerait en Physique par une « expérience de pensée » : on imagine que la procédure ou le processus -- entendu au sens commun du terme, pas nécessairement au sens probabiliste -- qui a généré les données peut être répété autant de fois qu'on le souhaite. Plus précisément, on postule que les données $X_1, X_2,\ldots,X_n$ -- pour simplifier, nous allons les supposer IID -- sont tirées d'une loi G de densité ou de distribution (suivant qu'elle est continue ou discrète) g. L'événement $E = (X_1 = x_1, X_2 = x_2,\ldots,X_n = x_n)$ se voit alors assigner une (densité de) probabilité : $$P(E) = P(X_1 = x_1, X_2 = x_2,\ldots,X_n = x_n) = P(x_1)P(x_2)\cdots P(x_n) = g(x_1) g(x_2) \cdots g(x_n)\, .$$ Dans le cadre que nous nous fixons, E n'est qu'un événement possible parmi un très grand nombre d'autres. Nous essayons ensuite d'expliquer les données $(x_1,x_2,\ldots,x_n)$ par une densité (distribution) que nous notons $f(x|θ)$ où $θ \in \mathbb{R}^p$. Si nous sommes de nature optimiste, nous pouvons même espérer qu'il existe un $θ_0$ tel que $f(x| θ_0) = g(x)$, mais nous n'avons pas besoin de le faire en général. La (densité de) probabilité de nos données devient alors: $$P(E|θ) = f(x_1|θ) f(x_2|θ)\cdots f(x_n|θ)\, .$$ Cela nous a amené à définir la vraisemblance par : $$L(θ) = c \cdot P(E|θ), \quad c > 0 \, .$$ La méthode d'estimation que nous avons privilégié est celle qui consiste à prendre : $$\hat{θ} = \arg \max_θ L(θ)\, ,$$ c'est la méthode du maximum de vraisemblance. La quantité $\hat{θ}$ étant l'estimateur du maximum de vraisemblance (EMV). Pour simplifier les calculs et améliorer nos évaluations numériques nous avons introduit la fonction de log-vraisemblance : $$l(θ) = \log L(θ)\, .$$ Pour nous guider dans l'interprétation des résultats, nous avons introduit la fonction de log-vraisemblance relative : $$r(θ) = l(θ) - l(\hat{θ}) \, .$$

Notre fonction de vraisemblance est explicitement définie comme une fonction des données (E), nous nous attendons donc à ce qu'elle « change » dans notre cadre théorique où E n'est qu'un événement possible parmi un très grand nombre d'autres. Il en va clairement de même pour $l(θ)$ et $r(θ)$. L'EMV $\hat{θ}$ étant définit comme l'argument auquel une fonction (aléatoire) atteind son maximum va, quant à lui, être une variable aléatoire. C'est-à-dire que dans notre cadre théorique, nous nous attendons à ce que $\hat{θ}$ change d'une répétition (tirage de E) à l'autre. Il devient alors naturel d'essayer de caractériser la loi de $\hat{θ}$ que nous allons désigner par loi d'échantillonage.

Exemple 11.1.1

Nous reprenons ici le premier exemple du chapitre 9 : 10 individus sont sélectionés aléatoirement à partir d'une grande population homogène et sont testés pour la tuberculose. Le but de ces mesures est d'estimer la proportion d'individus tuberculeux dans la population totale. Nous notons $x$ le nombre d'individus tuberculeux dans notre échantillon de taille 10. Nous avons vus que la log-vraisemblance peut alors s'écrire : $$l(θ) = x \log θ + (10-x) \log (1-θ), \quad \mathrm{pour} \quad 0 < θ < 1$$ et que l'EMV a pour expression : $$\hat{θ} = x/10, \quad \mathrm{pour} \quad 0 \le θ \le 1 \, .$$

Dans le cadre théorique de la section précédente, l'événement E se limite à $x$ et la loi $G$ est une loi binomiale, c'est-à-dire que nous avons : $$g(x) = f(x|θ_0) = \binom{10}{x} θ_0^x (1-θ_0)^{10-x} \, , \quad \mathrm{pour} \quad θ_0 \in [0,1] \, \quad \mathrm{et} \quad x = 0,1,\ldots,10 \, .$$ Pour chaque valeur possible de $θ_0$ nous avons donc 11 valeurs possibles de $x$ soient 11 fonctions de vraisemblance / log-vraisemblance / log-vraisemblance relative et 11 valeurs de $\hat{θ}$ possibles. Les 11 fonctions $r(θ)$ sont toujours les mêmes, ce n'est que leur fréquence relative qui va changer avec $θ_0$. Suivant Kalbfleisch, nous construisons les 11 fonctions $r(θ)$ possibles :

In [1]:
%pylab

Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib
In [2]:
from __future__ import print_function, division, unicode_literals, absolute_import

J'utilise maintenant une commande « magique » pour faciliter la génération d'un page HTML statique :

In [3]:
%matplotlib inline
In [29]:
theta_v = linspace(0.005,0.9995,201)
def fait_r_binom(x,n=10):
    theta_chapeau = x/n
    if x == 0 or x == 10:
       def r(theta): return x*log(theta)+(10-x)*log(1-theta)
    else:
        l_en_theta_chapeau = x*log(theta_chapeau)+(n-x)*log(1-theta_chapeau)
        def r(theta): return x*log(theta)+(n-x)*log(1-theta) - l_en_theta_chapeau
    return r
In [35]:
for x in range(11):
    r = fait_r_binom(x)
    plot(theta_v,r(theta_v),color='black',linestyle='solid')
    
axhline(log(0.1),color='red',linestyle='dashed')
xlabel(r'$\theta{}$')
ylabel(r'$r(\theta{})$')
ylim([-5,0])
Out[35]:
(-5, 0)

Suivant les conventions habituelles en statistique, Kalbfleisch introduit (page 97) le rapport de vraisemblance définit par : $$D = -2 r(θ_0) = -2 [l(θ_0) - l(\hat{θ})] \, .$$ Cette quantité n'a à première vue qu'un intérêt théorique puisqu'elle fait intervenir la « vraie » valeur du paramètre ($θ_0$) qui est précisément ce que nous cherchons. Nous pouvons néanmoins commencer à deviner son intérêt pratique en construisant sa distribution pour différentes valeurs de $θ_0$, par exemple pour 101 valeurs régulièrement espacées allant de 0 à 1. En effet, pour chaque valeur de $θ_0$ nous pouvons simplement obtenir la probabilité de chacune des 11 valeurs possibles de $x$ et donc celle de chacune des 11 valeurs possibles de l'EMV. Les graphes que nous venons de tracer nous donnent alors les valeurs possibles de $D$ avec leurs probabilités respectives. Nous commençons par définir une fonction qui calcule le rapport de vraisemblance pour $θ_0$, $x$ et $n$ fixés :

In [4]:
def D(theta_0,x,n):
    theta_chapeau = x/n
    if x == 0 or x == n: l_en_theta_chapeau = 0
    else: l_en_theta_chapeau = x*log(x)-x*log(n)+(n-x)*log(n-x)-(n-x)*log(n)
    return -2*(x*log(theta_0)+(n-x)*log(1-theta_0)-l_en_theta_chapeau)

Nous définissons ensuite une fonction qui nous renvoit un tuple à partir duquel la fonction de répartition du rapport de vraisemblance peut être obtenue pour $θ_0$ et $n$ fixés. Le tuple contient deux vecteurs de longueurs identiques ($n+1$) :

  • le premier contient les valeurs possibles et ordonnées (dans le sens croissant) du rapport de vraisemblance ;
  • le second contient les probabilité cumulées de ces valeurs.
In [5]:
def D_repartition(theta_0,n):
    x_v = range(n+1)
    p_x = array([math.factorial(n)/(math.factorial(x)*math.factorial(n-x))*theta_0**x*(1-theta_0)**(n-x) for x in x_v])
    D_v = array([D(theta_0,x,n) for x in x_v])
    idx = argsort(D_v)
    P = cumsum(p_x[idx])
    X = D_v[idx]
    return (X,P)

Ici nous employons deux fonctions du module numpy : argsort et cumsum ainsi que la fonction factorial du module math. En principe nous aurions pu écrire dans la deuxième ligne du corps de la fonction math.factorial(n)/math.factorial(x)/math.factorial(n-x) (c'est ce que j'avais fait au début), mais cela génère des problème quand $n$ devient grand : la division avec des flottants ne sait pas évaluer 1000!/50!.

Nous pouvons à présent calculer de façon systématique les fonctions de répartitions du rapport de vraisemblance pour différentes valeurs de $θ_0$ (de 0,1 à 0,5 par pas de 0,1 ; il est inutile d'aller au-delà de 0,5 par symétrie) et différentes tailles d'échantillons (10, 50 et 250) :

In [8]:
D_avec_10 = [D_repartition(theta,10) for theta in arange(1,6)/10]
D_avec_50 = [D_repartition(theta,50) for theta in arange(1,6)/10]
D_avec_250 = [D_repartition(theta,250) for theta in arange(1,6)/10]

figure(figsize=(12,8))
subplot(131)
step(D_avec_10[0][0],D_avec_10[0][1],where='post',color='red')
[step(Y[0],Y[1],where='post',color='black') for Y in D_avec_10[1:]]
xlim([0,10])
ylim([0,1])
xlabel(r'$D(\theta_0,x,n)$')
title('n=10')
ylabel(u'Probabilité')

subplot(132)
step(D_avec_50[0][0],D_avec_50[0][1],where='post',color='red')
[step(Y[0],Y[1],where='post',color='black') for Y in D_avec_50[1:]]
xlim([0,10])
ylim([0,1])
xlabel(r'$D(\theta_0,x,n)$')
title('n=50')

subplot(133)
step(D_avec_250[0][0],D_avec_250[0][1],where='post',color='red')
[step(Y[0],Y[1],where='post',color='black') for Y in D_avec_250[1:]]
xlim([0,10])
ylim([0,1])
xlabel(r'$D(\theta_0,x,n)$')
title('n=250')
Out[8]:
<matplotlib.text.Text at 0x7fe6abd15310>

Ici nous avons représenté en rouge la fonction de répartition correspondant à $θ_0 = 0,1$. Ce qui est frappant c'est que lorsque n croît, les fonctions de répartition deviennent toutes identiques, donc idépendantes du paramètre $θ_0$ de la loi qui a généré les données.

In [ ]: