Cours d'introduction à la statistique et à l'analyse de données : Notes de cours

De nombreuses « interfaces utilisateurs » sont disponible pour Python. La plus simple est démarrée en ligne de commande avec :

python

Pour pouvoir faire des calculs avec numpy et générer des graphiques avec matplotlib, le plus simple est alors d'importer dans leur intégralité quelques modules :

>>> from numpy import *
>>> from pylab import *

avant de passer en mode interactif avec :

>>> ion()

Ce qui permet de générer des graphes avec par exemple :

>>> plot([1,2,3,4,5],[5,2,7,0,1])

Pour éditer du code, vous devrez utiliser votre éditeur préféré (gEdit reconnaît le syntaxe Python, comme beaucoup d'autres éditeurs), puis copiez et coller dans votre console.

IDLE est l'environnement de développement intégré inclut par défaut dans toute distribution complète de Python (celle de Ubuntu n'est pas complète de ce point de vue et il faut installer IDLE en plus). IDLE est démarré comme Python avec :

idle

De même qu'avec l'interface de base, il faudra charger les modules numpy et pylab. Un différence est que les graphiques interactifs ne sont pas possible, il faut taper des séquences de commandes comme :

>>> figure(1)
>>> plot([1,2,3,4,5],[5,2,7,0,1])
>>> show()

et pendant que la fenêtre graphique est ouverte, il n'est pas possible d'entrer des commandes en console Python… Par contre, IDLE vient avec un éditeur assez performant de code Python.

Dans le cadre du cours, je vous recommande d'utiliser IPython, une console « hypertrophiée » pour faire du calcul scientifique interactif avec Python (qu'il faut installer en plus sur la plupart des distributions Linux). Comme les deux précédentes, on démarre avec :

ipython

Ensuite, la « commande magique » (la notion de commande magique fait partie du jargon, un peu prétentieux à mon sens, de IPython) :

%pylab

donne accès à numpy, matplotlib, etc.

Pour éditer du code, vous devrez utiliser votre éditeur préféré (gEdit reconnaît le syntaxe Python, comme beaucoup d'autres éditeurs), puis copiez et coller dans votre console avec la « commande magique » %paste.

L'éditeur emacs vient avec de nombreux modes Python. Personnellement, j'utilise le mode par défaut. Pour utiliser IPython à partir de emacs, il suffit d'affecter "ipython" à la variable python-shell-interpreter avec :

(setq python-shell-interpreter "ipython")

Pour avoir des choses un peu plus sophistiquées on pourra taper :

(setq
 python-shell-interpreter "ipython2"
 python-shell-prompt-regexp "In \\[[0-9]+\\]: "
 python-shell-prompt-output-regexp "Out\\[[0-9]+\\]: "
 python-shell-completion-setup-code
   "from IPython.core.completerlib import module_completion"
 python-shell-completion-module-string-code
   "';'.join(module_completion('''%s'''))\n"
 python-shell-completion-string-code
   "';'.join(get_ipython().Completer.all_completions('''%s'''))\n")

Ce qui précède est bien sûr loin d'être exhaustif, pour plus de choix, consultez le wiki Python. Si vous cherchez des interfaces sophistiquées type environnement de développement intégré (Integrated Development Environment IDE), jetez un œil sur pycharm — dont la version « de base » est gratuite — et sur spyder.

Nous allons considérer la densité suivante : \[f(x) = \frac{1}{2}\, \phi(x;0,1) + \frac{1}{10} \, \sum_{j=0}^4 \phi\left(x;(j/2)-1,1/10\right) \, ,\] où \(\phi(x; \mu, \sigma)\) désigne la densité d'une loi normale de moyenne \(\mu\) et d'écart type \(\sigma\).

import scipy
from scipy.stats import norm
def bart_simpson_pdf(x):
    res = 0.5*norm.pdf(x)
    for j in range(5):
        res += norm.pdf(x,loc=j*0.5-1,scale=0.1)*0.1
    return res

def bart_simpson_rvs(size=1):
    weights = array([0.5,0.1,0.1,0.1,0.1,0.1])
    means = array([0,-1,-0.5,0,0.5,1])
    sigmas = array([1,0.1,0.1,0.1,0.1,0.1])
    res = norm.rvs(size=size)
    k = list(numpy.random.choice(arange(6),size=size,p=weights))
    locs = means[k]
    scales = sigmas[k]
    res *= scales
    res += locs
    return res

seed(20110928)
bart_simpson_test0 = bart_simpson_rvs(1000)

def bart_simpson_cdf(x):
    res = 0.5*norm.cdf(x)
    for j in range(5):
        res += norm.cdf(x,loc=j*0.5-1,scale=0.1)*0.1
    return res

def fait_fre(echantillon):
    y = sort(echantillon)
    n = float(len(echantillon))
    def fre(x):
        return sum(y <= x)/n
    return fre

bart_simpson_test0_fre = fait_fre(bart_simpson_test0)
bb = arange(-3,3.0001,0.0001)
bart_simpson_test0_delta = array([bart_simpson_test0_fre(bb[i+1])-bart_simpson_test0_fre(bb[i]) for i in range(len(bb)-1)])
bb_mid = bb[:-1] + diff(bb)/2.
kernelBest = norm.pdf(arange(-0.25,0.2501,0.0001),scale=0.05)
from scipy.signal import fftconvolve
bart_simpson_test0_pdf_best = fftconvolve(bart_simpson_test0_delta,kernelBest,mode='same')
kernelBest_over_10 = norm.pdf(arange(-0.025,0.02501,0.0001),scale=0.005)
bart_simpson_test0_pdf_best_over_10 = fftconvolve(bart_simpson_test0_delta,kernelBest_over_10,mode='same')
kernelBest_times_10 = norm.pdf(arange(-2.5,2.501,0.0001),scale=0.5)
bart_simpson_test0_pdf_best_times_10 = fftconvolve(bart_simpson_test0_delta,kernelBest_times_10,mode='same')

La figure suivante montre la densité de Bart Simpson ainsi que trois estimations de celle-ci basée sur un estimateur à noyau construit à partir du même échantillon de taille 1000 en utilisant trois valeurs du paramètre de lissage:

subplot(221)
xx = linspace(-3.,3,501)
plot(xx,bart_simpson_pdf(xx),color='black')
title(u'Vraie densité')
xlim([-3,3])
ylim([0,1])
subplot(222)
plot(bb_mid,bart_simpson_test0_pdf_best_over_10,color='black')
title(u'Pas assez lissée')
xlim([-3,3])
ylim([0,1])
subplot(223)
plot(bb_mid,bart_simpson_test0_pdf_best,color='black')
title(u"C'est bon !")
xlim([-3,3])
ylim([0,1])
subplot(224)
plot(bb_mid,bart_simpson_test0_pdf_best_times_10,color='black')
title(u'Trop lissée')
xlim([-3,3])
ylim([0,1])
savefig('img/coursIntroStats2014/Clone-Fig6p1-AONPS.png')
close()
'img/coursIntroStats2014/Clone-Fig6p1-AONPS.png'

Nous allons juger la qualité d'un estimateur \(\widehat{f_n}\) par son erreur quadratique intégréee moyenne, \(R = \mathrm{E}(L)\) où : \[L = \int \left(\widehat{f_n}(x) - f(x)\right)^2 dx\] est l'erreur quadratique intégrée. Ici l'espérance est prise par rapport à la loi d'un vecteur de \(n\) variables aléatoires IID de loi \(\mathrm{F}\) ; nous notons \(\widehat{f_n}(x)\) ce que nous pourrions écrire plus explicitement : \(\widehat{f_n}(x;X_1,\ldots,X_n)\). Ainsi avec des notations plus lourdes nous pouvons écrire : \[R = \int_{X_1} \ldots \int_{X_n} \left(\int \left(\widehat{f_n}(x;X_1,\ldots,X_n) - f(x)\right)^2 dx\right) dF(X_1) \cdots dF(X_n) \; .\]

La méthode habituelle pour estimer l'erreur quadratique intégrée moyenne est la validation croisée ou, plus précisemment, ce qu'on désigne en anglais par leave-one-out-cross-validation. Avant de la définir, nous allons un peu « travailler » notre expression de l'erreur quadratique intégrée : \[\begin{array}{l,l,l}L(h) & = & \int \left(\widehat{f_n}(x) - f(x)\right)^2 dx \, , \\ & = & \int \widehat{f_n}^2(x)\, dx - 2\, \int \widehat{f_n}(x) \, f(x) dx + \int f^2(x) dx \, .\end{array}\] Comme le dernier terme ne dépend pas de \(h\), minimiser l'erreur quadratique intégrée moyenne revient à minimiser l'espérance de : \[J(h) = \int \widehat{f_n}^2(x)\, dx - 2\, \int \widehat{f_n}(x) \, f(x) dx \, .\] Pour faire cours, nous allons appeler \(J(h)\) l'erreur quadratique intégrée moyenne, même si cette quantité diffère de la vraie erreur quadratique intégrée moyenne par le terme constant \(\int f^2(x) dx\).

Notre estimateur de l'erreur quadratique intégrée moyenne par validation croisée est alors définit par : \[\widehat{J}(h) = \int \widehat{f_n}^2(x)\, dx - \frac{2}{n} \sum_{i=1}^n \widehat{f}_{(-i)}(X_i) \] où \(\widehat{f}_{(-i)}\) est l'estimateur de densité obtenu après avoir retiré la \(i\) -ième observation de l'échantillon. Nous appellerons \(\widehat{J}(h)\) le « score de validation croisée » ou « l'erreur quadratique intégrée moyenne estimée ».

L'estimateur de densité le plus simple et, sans aucun doute, le plus communément employé est l'histogramme. Supposons que \(f\) admet pour support un intervalle que nous pourrons prendre (moyennant une transformation linéaire si nécessaire) comme étant \([0,1]\). Soit \(m\) un entier, nous définissons les classes (bins en anglais) : \[C_1 = \left[0,\frac{1}{m}\right), \quad C_2 = \left[\frac{1}{m},\frac{2}{m}\right), \quad \ldots, \quad C_m = \left[\frac{m-1}{m},1\right) \, .\] Nous définissons :

• la largeur de classe (binwidth) \(h = 1/m\) ;
• \(Y_j\) la variable aléatoire décrivant le nombre d'observations dans la classe \(C_j\) ;
• \(\hat{p}_j = Y_j/n\) ;
• \(p_j = \int_{C_j} \, f(u) \, du\).

L'(estimateur) histogramme est alors définit par : \[\widehat{f_n}(x) = \sum_{j=1}^m \frac{\hat{p}_j}{h}\, I(x \in C_j) \, .\]

Pour comprendre l'origine de cet estimateur, il suffit de remarquer que : \[\mathrm{E}\left(\widehat{f_n}(x)\right) = \frac{\mathrm{E}(\hat{p}_j)}{h} = \frac{p_j}{h} = \frac{\int_{C_j} \, f(u) \, du}{h} \approx \frac{f(x)\, h}{h} = f(x) \, .\]

On en déduit alors :

Théorème Soient \(x\) et \(m\) fixés et soit \(C_j\) la classe contenant \(x\) ; alors \[\mathrm{E}\left(\widehat{f_n}(x)\right) = \frac{p_j}{h} \quad \mathrm{et} \quad \mathrm{Var}\left(\widehat{f_n}(x)\right) = \frac{p_j\, (1-p_j)}{n\, h^2} \, .\]

Le théorème suivant montre qu'il est inutile de (re)calculer \(n\) estimateurs pour obtenir la valeur de l'erreur quadratique intégrée moyenne estimée.

Théorème Nous avons l'identité suivante : \[\widehat{J}(h) = \frac{2}{h\, (n-1)} - \frac{n+1}{h\, (n-1)} \sum_{j=1}^m \hat{p}_j^2 \, .\]

Preuve (Avant de lire les lignes suivantes essayer d'obtenir la preuve vous même.) Par définition : \[\widehat{J}(h) = \int \widehat{f_n}^2(x)\, dx - \frac{2}{n} \sum_{i=1}^n \widehat{f}_{(-i)}(X_i) \, , \] \[\widehat{f_n}(x) = \sum_{j=1}^m \frac{\hat{p}_j}{h}\, I(x \in C_j) \, ,\] \[\hat{p}_j = Y_j/n \, .\] Le premier terme de \(\widehat{J}(h)\) devient donc : \[\begin{array}{l,l,l} \int \widehat{f_n}^2(x)\, dx & = & \int \left(\sum_{j=1}^m \frac{\hat{p}_j}{h}\, I(x \in C_j)\right) \, \left(\sum_{k=1}^m \frac{\hat{p}_k}{h}\, I(x \in C_k)\right) \, dx \\ & = & \sum_{j=1}^m \sum_{k=1}^m \frac{\hat{p}_j}{h} \frac{\hat{p}_k}{h} \int I(x \in C_j) \, I(x \in C_k) \, dx \\ & = & \sum_{j=1}^m \frac{\hat{p}^2_j}{h} \end{array}\] où on a utlisé le fait : \[I(x \in C_j) \, I(x \in C_k) = \left\{ \begin{array}{l l} I(x \in C_j) \quad \mathrm{si} \quad j=k \, ,\\ 0 \quad \mathrm{sinon} \end{array}\right.\] et : \[\int I(x \in C_j)\, dx = \int_{\frac{j-1}{m}}^{\frac{j}{m}} dx = \frac{1}{m} = h \, .\]

Pour calculer le second terme, \(\sum_{i=1}^n \widehat{f}_{(-i)}(X_i)\), nous commençons par trouver l'expression de \(\widehat{f}_{(-i)}(x)\) : \[\begin{array}{l,l,l} \widehat{f}_{(-i)}(x) & = & \sum_{j=1}^m \frac{\hat{p}_{(-i),j}}{h}\, I(x \in C_j) \\ & = & \sum_{j=1}^m \frac{Y_{(-i),j}}{h \, (n-1)}\, I(x \in C_j) \\ & = & \sum_{j=1,j \neq k}^m \frac{Y_{(-i),j}}{h \, (n-1)}\, I(x \in C_j) + \frac{Y_{(-i),k}}{h \, (n-1)}\, I(x \in C_k) \end{array}\]

où \(k\) est l'indice de la classe contenant \(X_i\). Mais : \[Y_{(-i),j} = \left\{ \begin{array}{l l} Y_j \quad \mathrm{si} \quad j \neq k \\ Y_j - 1 \quad \mathrm{si} \quad j = k\end{array} \right.\] d'où : \[\begin{array}{l,l,l} \widehat{f}_{(-i)}(x) & = & \frac{n}{n-1} \sum_{j=1}^m \frac{Y_j}{h \, n}\, I(x \in C_j) - \frac{1}{h\, (n-1)} I(x \in C_k) \\ & = & \frac{n}{n-1} \sum_{j=1}^m \frac{\hat{p}_j}{h}\, I(x \in C_j) - \frac{1}{h\, (n-1)} I(x \in C_k) \end{array}\]

On en déduit que : \[\widehat{f}_{(-i)}(X_i) = \frac{n}{n-1} \frac{\hat{p}_k}{h} - \frac{1}{h\, (n-1)}\] et que : \[\begin{array}{l,l,l} \sum_{i=1}^n \widehat{f}_{(-i)}(X_i) & = & \frac{n}{n-1} \sum_{k=1}^m\frac{Y_k \, \hat{p}_k}{h} - \frac{n}{h\, (n-1)} \\ & = & \frac{n^2}{n-1} \sum_{k=1}^m\frac{\hat{p}^2_k}{h} - \frac{n}{h\, (n-1)} \end{array}\] où le fait que \(Y_k\) des \(X_i\) se trouvent dans la classe \(k\) a été utilisé. En sommant les deux termes en en réarrangeant il vient : \[\widehat{J}(h) = \frac{2}{h\, (n-1)} - \frac{n+1}{h\, (n-1)} \sum_{j=1}^m \hat{p}_j^2 \, .\]

Question Pourquoi retire-t-on l'observation de \(X_i\) dans notre définition de l'erreur quadratique intégrée moyenne estimée ?

Réponse Pour le comprendre, vous pourrez considérer l'estimateur : \[\check{J}(h) = \int \widehat{f_n}^2(x)\, dx - \frac{2}{n} \sum_{i=1}^n \widehat{f}_{n}(X_i) \, .\] Avec des \(X_i\) tous différents, il est possible de trouver un \(h < \min \left(X_{(i+1)}-X_{(i)}\right)\), c'est-à-dire que pour \(h\) assez petit, les classes de notre histogrammes contiennent soit 0 soit 1 observation (\(Y_i = 0\) ou \(Y_i = 1\) pour tout \(i\)). On montre alors facilement que pour tout \(h\) assez petit : \[\check{J}(h) = - \frac{1}{h \, n} \, .\] C'est-à-dire que le \(h\) optimal est nul et l'histogramme devient un peigne de Dirac : \[\check{f}_n(x) = \frac{1}{n} \sum_{i=1}^n \delta_{X_i}(x) \, .\] En conséquence, si \(X_{extra}\) est une nouvelle observation de \(F\) nous aurons : \(\check{f}_n(X_{extra}) = 0\)…

Supposons que les v.a. \(X_1,\ldots,X_n\) sont IID de loi uniforme continue sur l'intervalle \([0,\theta]\) pour un \(\theta > 0\). La fonction de vraisemblance est : \[\mathcal{L}(\theta) = \theta^{-n} I(0 \le x_1,\ldots,x_n \le \theta) = \theta^{-n} I\left(\theta \ge \max(x_1,\ldots,x_n)\right)\, .\] Donc \(\mathcal{L}(\theta) = 0\) si \(\theta < \max(x_1,\ldots,x_n)\) alors que \(\mathcal{L}(\theta)\) est strictement décroissante si \(\theta \ge \max(x_1,\ldots,x_n)\). Donc \(\mathcal{L}(\theta)\) est maximale en \(\theta = \max(x_1,\ldots,x_n)\) et : \[\hat{\theta} = X_{(n)} = \max(x_1,\ldots,x_n) \, .\]

Supposons que les v.a. \(X_1,\ldots,X_n\) sont IID de loi de Poisson de paramètre \(\lambda > 0\). La fonction de vraisemblance est : \[\mathcal{L}(\lambda) = \prod_{i=1}^n \frac{\lambda^{x_i} \, \exp(-\lambda)}{x_i!}\] et la log-vraisemblance est : \[\ln \mathcal{L}(\lambda) = -n \, \lambda + \ln(\lambda) \, \sum_{i=1}^n x_i - \sum_{i=1}^n \ln(x_i!) \, .\] En supposant \(\sum_{i=1}^n x_i > 0\) et en dérivant par-rapport à \(\lambda\) il vient : \[\frac{d}{d\lambda} \ln \mathcal{L}(\lambda) = -n + \frac{1}{\lambda} \sum_{i=1}^n x_i\] et la dérivée s'annule en \(\lambda = \overline{x}\). Nous vérifions maintenant que ce point correspond bien à un maximum en calculant la dérivée seconde : \[\frac{d^2}{d\lambda^2} \ln \mathcal{L}(\lambda) = - \frac{1}{\lambda^2} \sum_{i=1}^n x_i \, ,\] qui est toujours négative. La fonction de vraisemblance atteint donc son unique maximum en \(\overline{x}\) et l'EMV de \(\lambda\) est \(\hat{\lambda} = \overline{X}\), dès lors que \(\sum_{i=1}^n X_i > 0\). Si \(\sum_{i=1}^n X_i = 0\), l'EMV n'existe pas.

La fonction génératrice des moments d'une variable aléatoire \(X\) est définie par : \[M_X(t) = \mathrm{E}\left(\exp(tX)\right),\quad |t| \le b > 0\, ,\] lorsque \(M_X(t) < \infty\).

Comme son nom l'indique, la fonction génératrice des moments est employée pour obtenir les moments de la v.a. \(X\) : \(\mathrm{E}\left(X\right), \mathrm{E}\left(X^2\right), \mathrm{E}\left(X^3\right), \ldots\). \(\mathrm{E}\left(X^k\right)\) est, par définition, le $k$-ième moment de \(X\). Nous avons alors le théorème suivant :

Théorème Soit \(X\) une variable aléatoire de fonction génératrice des moments \(M_X(t)\). Alors pour tout \(r > 0\), \(\mathrm{E}\left(|X|^r\right) < \infty\) et \(\mathrm{E}\left(X\right) = M'_X(0)\), \(\mathrm{E}\left(X^2\right) = M''_X(0)\) et plus généralement, \(\mathrm{E}\left(X^k\right) = M^{(k)}_X(0)\), où \(M^{(k)}_X\) désigne la $k$-ième dérivée de \(M_X\).

Supposons que les v.a. \(X_1,\ldots,X_n\) sont IID de loi géométrique dont la fonction de masse est : \[f(x;\theta) = \theta \, (1-\theta)^{x-1}, \quad x=1,2,\ldots, \quad 0 < \theta < 1\; .\] Attention lors de la présentation du cours, j'ai suivi le bouquin de Keith Knight et utilisé une définition moins habituelle de la fonction de masse : \(f(x;\theta) = \theta \, (1-\theta)^x, \quad x=0,1,2,\ldots, \quad 0 < \theta < 1\). La suite est en fait plus simple avec la définition plus « classique ».

Vous pouvez « voir » à quoi cette loi correspond en considérant que \(\theta\) est la probabilité de gagner au tirage d'une loterie dont les tirages successifs sont indépendants et identiquement distribués, \(f(x;\theta)\) est alors la probabilité de gagner pour la première fois au $x$-ième tirage. Le nom de la loi vient de la série géométrique de raison \(q\), pour laquelle, si vous vous souvenez bien : \[\sum_{i=0}^{\infty} q^i = \frac{1}{1-q}, \quad |q| < 1 \; .\] En posant \(q = (1-\theta)\) vous voyez immédiatement que : \[\sum_{x=1}^{\infty} f(x;\theta) = 1\; .\]

Ici la vraisemblance de \(\theta\) étant donné \(X_1=x_1,\ldots,X_n=x_n\) s'écrie : \[\mathcal{L}(\theta) = \prod_{i=1}^n \theta \, (1-\theta)^{x_i-1}\] et la log-vraisemblance est : \[\mathcal{l}(\theta) = \log \mathcal{L}(\theta) = n\, \log(\theta) + \log(1-\theta) \, \sum_{i=1}^n (x_i-1) \; .\] En dérivant deux fois on établit que l'EMV de \(\theta\) est : \[\hat{\theta} = \frac{1}{\overline{X}_n} \; ,\] où \[\overline{X}_n \doteq \frac{1}{n} \sum_{i=1}^n X_i \; .\] L'EMV est ainsi une fonction continue de la moyenne empirique et nous avons : \[\overline{X}_n \stackrel{P}{\rightarrow} \mathrm{E}(X) \quad \textrm{(loi des grands nombres)} \, .\] Comme \(\hat{\theta} = g(\overline{X}_n)\) où \(g(x)=1/x\) est une fonction continue (pour \(x>0\), ce qui est le cas ici), la convergence d'une fonction d'une variable aléatoire nous donne : \[\hat{\theta} \stackrel{P}{\rightarrow} \frac{1}{\mathrm{E}(X)}\] ce qui prouve la convergence de l'EMV.

Le théorème central limite (TCL) nous permet d'affirmer que : \[\sqrt{n} \, \left(\overline{X}_n - \mathrm{E}(X)\right) \stackrel{L}{\rightarrow} \mathcal{N}(0,\sigma^2_X)\] si \(\sigma^2_X < \infty\) et la méthode delta que si \(g\) est dérivable telle que \(g'\left(\mathrm{E}(X)\right) \neq 0\) alors : \[\sqrt{n} \, \left(g\left(\overline{X}_n\right) - g\left(\mathrm{E}(X)\right)\right) \stackrel{L}{\rightarrow} \mathcal{N}\left[0,\sigma^2_X g'\left(\mathrm{E}(X)\right)^2\right] \, .\]

On pourra préférer écrire ces deux dernières équations : \[\overline{X}_n \stackrel{L}{\rightarrow} \mathrm{E}(X) + \frac{\sigma_X}{\sqrt{n}} \, Z\; ,\] où \(Z\) désigne une v.a. de loi normale centrée réduite et : \[g\left(\overline{X}_n\right) \stackrel{L}{\rightarrow} g\left(\mathrm{E}(X)\right) + |g'\left(\mathrm{E}(X)\right)|\,\frac{\sigma_X}{\sqrt{n}} \, Z\; .\]

Il nous reste donc à trouver \(\mathrm{E}(X)\) et \(\sigma_X\) pour \(X\) v.a. de loi géométrique de paramètre \(\theta\). Pour ce faire, nous pouvons, comme nous l'avons vu en cours, employer une méthode directe qui repose essentiellement sur l'identité : \[S(\theta) = \sum_{x=1}^{\infty} \theta \, (1-\theta)^{x-1} = 1 \, ,\] suivie d'une dériavtion par-rapport à \(\theta\) en interchangeant la dérivation avec la sommation, puis, après quelques manipulations, une identification. Clairement, cette méthode, si elle peut toujours être mise en œuvre, prend du temps et nous expose aux erreurs de calculs. Une autre solution est d'aller voir le site wikipédia de la loi géométrique où on trouvera : \[\mathrm{E}(X) = 1/\theta \textrm{ et } \sigma_X = \sqrt{1-\theta}/\theta\, .\] Faites attention ici à la convention adoptée lors de la définition de la loi (la « classique » avec l'indice qui commence à 1 ou celle, plus rare, de Knight dont l'indice commence à 0). Les « bons bouquins », comme ceux de Rice et Wasserman, ont une table des lois usuelles avec leurs propriétés comme la moyenne, l'écart type, etc… On peut aussi faire appel au module SymPy de Python, module dédié à la résolution symbolique de problèmes mathématiques comme l'intégration, la différentiation, le calcul des limites, etc… Nous allons l'employer ici pour faire des calculs symboliques à partir de la fonction génératrice des moments, cette dernière est ici : \[M_X(t) = \mathrm{E}\left(\exp(t X)\right) = \sum_{x=1}^{\infty} \exp(t x) \, \theta \, (1-\theta)^{x-1}\; ,\] ce qui donne quasi immédiatement : \[M_X(t) = \frac{\theta}{1-\theta} \sum_{x=1}^{\infty} \left(\exp(t) \, (1-\theta)\right)^x = \frac{\theta \, \exp t}{1 - \exp(t) \, (1-\theta)}\, \textrm{ pour } t < -\log(1-\theta)\; .\] Nous définissons cette fonction génératrice des moments de la façon suivante — avec SymPy nous devons introduire des « variables symboliques » qui restent non-évaluées :

from sympy import *
theta = Symbol('theta', real=True, positive=True)
t = Symbol('t')
M_X = theta * exp(t)/(1-exp(t)*(1-theta))

Si vous souhaitez ne pas importer tout l'espace de noms comme nous l'avons fait ici (avec from symport *) vous pouvez entrer quelque chose comme :

import sympy as sy
theta = sy.Symbol('theta', real=True, positive=True)
t = sy.Symbol('t')
M_X = theta * sy.exp(t)/(1-sy.exp(t)*(1-theta))

Remarquez que nous appelons alors la version « symbolique » de l'exponentielle (sy.exp) et non la version numérique qui est appelée par défaut (exp). Il en va alors de même pour toutes les fonctions classiques qui constituent des expressions symboliques.

Pour obtenir la dérivée première nous utilisons la fonction diff (ou sy.diff si nous n'avons pas importé tout l'espace de noms) :

diff(M_X,t)

theta*(-theta + 1)*exp(2*t)/(-(-theta + 1)*exp(t) + 1)**2 + theta*exp(t)/(-(-theta + 1)*exp(t) + 1)

Quand les résultats sont « un peu longs », leur appliquer la fonction simplify est en général une bonne idée :

simplify(diff(M_X,t))

theta*exp(t)/((theta - 1)*exp(t) + 1)**2

Maintenant, pour obtenir la valeur de notre premier moment, nous devons remplacer t par 0 ou substituer 0 à t en langage de calcul symbolique. Cette opération est effectuée en appelant la méthode subs sur l'objet précédent, de la façon suivante :

simplify(diff(M_X,t)).subs(t,0)

1/theta

Ces opérations de substitutions constituent la base du calcul symbolique et vous êtes invités à lire le paragraphe du tutoriel qui s'y rapporte.

Remarque : les sorties les plus simples générées par SymPy sont montrées ici (de l'ASCII pur). Vous pouvez générer des sorties plus sophistiquées avec la commande :

init_printing()

Pour une description plus complète des différentes sorties (et sur comment les obtenir), consultez la page de la documentation correspondante.

La variance est le second moment auquel on soustrait le carré du premier (la moyenne), ce qu'on obtient avec :

simplify((diff(M_X,t,2)-diff(M_X,t,1)**2).subs(t,0))

(-theta + 1)/theta**2

En fin de compte, nous avions établi que : \[\hat{\theta} \stackrel{P}{\rightarrow} \frac{1}{\mathrm{E}(X)}\; ,\] or nous venons de voir que : \(\mathrm{E}(X) = 1/\theta\), ce qui nous donne : \[\hat{\theta} \stackrel{P}{\rightarrow} \theta \; .\] La loi normale asymptotique : \[\sqrt{n} \, \left(g\left(\overline{X}_n\right) - g\left(\mathrm{E}(X)\right)\right) \stackrel{L}{\rightarrow} \mathcal{N}\left[0,\sigma^2_X g'\left(\mathrm{E}(X)\right)^2\right] \, ,\] devient, en utilisant le fait que \(g(x) = 1/x\) et \(g'(x) = -1/x^2\) (attention il faut prendre ici \(x = 1/\theta\)): \[\sqrt{n} \, \left(\hat{\theta} - \theta\right) \stackrel{L}{\rightarrow} \mathcal{N}\left[0,(1-\theta) \, \theta^2\right] \, .\]

Il peut être intéressant à ce stade de comparer l'approximation asymptotique ci-dessus avec la distribution exacte puisque cette dernière est disponible. En effet, notons \(T = \sum_{i=1}^s X_i\), où les \(X_i\) sont des v.a. IID de loi géométrique de paramètre \(\theta\). Nous avons clairement : \(\overline{X}_s = T/s\). Pour que \(T=t\) il faut que parmi les \(t-1\) tirages qui précèdent, nous aillons \(s-1\) « succès ». Rappelez-vous que chaque tirage suivant une loi géométrique de paramètre \(\theta\) peut être vu comme une suite de \(x\) tirages IID suivant une loi de Bernoulli pour laquelle \(\theta\) est la probabilité de succès et où \(x=1,2,\ldots\) est la tirage pour lequel le premier succès est observé. Ainsi, \(s\) tirages indépendants suivant une loi géométriques peuvent être vus comme une suite de tirages de suivant une loi de Bernoulli jusqu'à ce que le $s$-ième succès soit observé. On obtient ainsi la fonction de masse de la v.a. \(T\) : \[\mathrm{Pr}\{T=t\} = {t-1 \choose s-1}\theta^s \, (1-\theta)^{t-s} \; , \quad t = s,s+1,\ldots\] La v.a. \(T\) suit une loi binomiale négative de paramètres \(s=1,2,\ldots\) et \(0<\theta<1\) (on peut en fait généraliser la définition à tout \(n>0\) réel en utilisant la fonction gamma à la place des factoriels dans l'équation précédente). Attention : la loi est parfois définie (comme sur le site wikipédia, où vous trouverez aussi une explication du nom de la loi) comme donnant la probabilité du nombre d'échecs observés jusqu'au $s$-ième succès, on a alors : \[\mathrm{Pr}\{E=e\} = {s+e-1 \choose s-1}\theta^s \, (1-\theta)^{e} \; , \quad e=0,1,\ldots \; ,\] les paramètres de la loi sont les mêmes qu'avec la définition précédente.

La convergence en loi que nous pouvons réécrire : \[\sqrt{\frac{n}{(1-\theta) \, \theta^2}} \, \left(\hat{\theta} - \theta\right) \stackrel{L}{\rightarrow} \mathcal{N}\left(0,1\right) \, ,\] nous dit que la fonction de répartition de la v.a. : \[\hat{\eta}_n \doteq \sqrt{\frac{n}{(1-\theta) \, \theta^2}} \, \left(\hat{\theta} - \theta\right)\; ,\] (je rajoute l'indice pour garder trace de la taille de l'échantillon) admet pour limite quand \(n\) tend vers l'infini, la fonction de répartition d'une loi normale centrée réduite. Mais \(\hat{\theta}_n = n/T_n\) donc : \[\begin{array}{l,l,l} F_{\hat{\theta}_n}(x) & = & \mathrm{Pr}\{\hat{\theta}_n \le x\} \textrm{ définition de la fonction de répartition,}\\ & = & \mathrm{Pr}\{n/T_n \le x\} \textrm{ définition de l'EMV,} \\ & = & \mathrm{Pr}\{T_n \ge n/x\} \, .\end{array}\] Et \(T_n\) le nombre total de tirages de Bernoulli jusqu'à l'obtention du $n$-ième succès peux aussi s'écrire : \(T_n = E_n + n\), où \(E_n\) est le nombre d'échecs jusqu'au $n$-ième succès, d'où : \[\begin{array}{l,l,l} F_{\hat{\theta}_n}(x) & = & \mathrm{Pr}\{T_n \ge n/x\} \\ & = & \mathrm{Pr}\{E_n + n \ge n/x\} \\ & = & \mathrm{Pr}\{E_n \ge n/x - n\} \\ & = & 1 - \mathrm{Pr}\{E_n \le n/x - n\} \\ & = & 1 - F_{E_n}\left(n/x - n\right)\, .\end{array}\] J'ai fait toute cette « gymnastique » par-ce que le site wikipédia me dit que la fonction de répartition d'une loi binomiale négative est une fonction bêta incomplète régularisée, c'est-à-dire que : \[F_{E_n}(y;n,\theta) = I_{\theta}(n,y+1) = \frac{\int_0^{\theta}t^{n-1} (1-t)^y dt}{\int_0^{1}t^{n-1} (1-t)^y dt} \; .\] « Heureusement » pour nous, le sous module special de scipy contient cette fonction et nous pouvons définir ainsi une fonction Python qui renvoie la valeur de la fonction de répartition de notre \(\hat{\theta}_n\) de la façon suivante (de ce qui suit, comme d'habitude, je suppose que vous utilisez IPython et que la commande magique %pylab a été exécutée) :

from scipy.special import betainc
def F_theta_hat(x,n,theta): 
    return 1-betainc(n,n/x-n,theta)

Nous obtenons maintenant facilement, la fonction de répartition de \(\hat{\eta}_n\) (qui doit être asymptotiquement normale, centrée, réduite) à partir de celle de \(\hat{\theta}\) et de la définition de \(\hat{\eta}_n\) : \[\begin{array}{l,l,l} F_{\hat{\eta}_n}(x) & = & \mathrm{Pr}\{\hat{\eta}_n \le x\} \\ & = & \mathrm{Pr}\left\{\sqrt{\frac{n}{(1-\theta) \, \theta^2}} \, \left(\hat{\theta}_n - \theta\right) \le n/x \right\} \\ & = & \mathrm{Pr}\left\{\hat{\theta}_n \le \theta \, \sqrt{\frac{1-\theta}{n}}\, x + \theta \right\} \\ & = & F_{\hat{\theta}_n}\left(\theta \, \sqrt{\frac{1-\theta}{n}}\, x + \theta\right)\, .\end{array}\]

D'où la définition de notre fonction de répartition en Python :

def F_eta_hat(x,n,theta):
    return F_theta_hat(theta*(sqrt((1-theta)/n)*x+1),n,theta)

Théorème Soient \(X_1,\ldots,X_n\) des v.a. IID de densité (ou « fonction de masse ») \(f_0 = f(x;\theta_0) \in \mathcal{F}\) où \(\mathcal{F} = \{f(x;\theta) : \theta \in \Theta\}\) est notre ensemble de densités « candidates ». Définissons : \[M_n(\theta) \doteq \frac{1}{n} \sum_{i=1}^{n} \ln \frac{f(X_i;\theta)}{f(X_i;\theta_0)}\] et \[M(\theta) \doteq \mathrm{E}_0\{\ln \frac{f(X_i;\theta)}{f(X_i;\theta_0)}\} = -D(\theta_0,\theta) \; .\] Si \[\sup_{\theta \in \Theta}|M_n(\theta) - M(\theta)| \stackrel{P}{\rightarrow} 0\] et si pour tout \(\epsilon > 0\), \[\sup_{\theta : |\theta - \theta_0 | \ge \epsilon} M(\theta) < M(\theta_0) \; ,\] alors : \(\hat{\theta} \stackrel{P}{\rightarrow} \theta_0\), où \(\hat{\theta}\) désigne l'EMV.

Preuve Par définition de \(\hat{\theta}\), \(M_n(\theta)\) est maximale en \(\hat{\theta}_n\) donc \(M_n(\hat{\theta}_n) \ge M_n(\theta_0)\). D'où \[\begin{array}{l,l,l} M(\theta_0) - M(\hat{\theta}_n) & = & M_n(\theta_0) - M(\hat{\theta}_n) + M(\theta_0) - M_n(\theta_0) \, , \\ & \le & M_n(\hat{\theta}_n) - M(\hat{\theta}_n) + M(\theta_0) - M_n(\theta_0) \, , \\ & \le & \sup_{\theta \in \Theta} |M_n(\theta) - M(\theta)| + M(\theta_0) - M_n(\theta_0) \, , \\ & \stackrel{P}{\rightarrow} & 0 \, .\end{array}\] Ainsi, quelque soit \(\delta > 0\) : \[\mathrm{Pr}\left\{M(\hat{\theta}_n) < M(\theta_0) - \delta \right\} \rightarrow 0 \; .\] Choisissons un \(\epsilon > 0\) quelconque. Notre dernière hypothèse nous dit qu'il existe un \(\delta > 0\), tel que \(|\theta - \theta_0| \ge \epsilon\) implique \(M(\theta) < M(\theta_0) - \delta\). Donc \[\mathrm{Pr}\left\{|\hat{\theta}_n - \theta_0| \ge \epsilon\right\} \le \mathrm{Pr}\left\{M(\hat{\theta}_n) < M(\theta_0) - \delta \right\} \rightarrow 0 \; .\] Ce qui complète la preuve.

Remarques :

• lors du cours, nous avons vu une condition un peu moins générale que la seconde hypothèse ci-dessus (qui revient à supposer l'existence d'un maximum global à l'intérieur de \(\Theta\)) pour la convergence puisque, suivant le bouquin de Keith Knight, nous avions supposé que la suite de fonctions \(M_n\) était une suite de fonctions concaves (j'ai suivi ici le bouquin de Wasserman) ;
• pour pouvoir appliquer le théorème, il faut vérifier une condition stochastique (convergence uniforme des \(M_n\) vers \(M\)), mais cette condition est relativement facile à vérifier dès lors que les données sont IID, il faut vérifier ensuite une condition déterministe, qui dans le cas d'un paramètre scalaire peut se faire en construisant le graphe de la fonction.

Pour \(f(X;\theta)\) une densité de probabilité (ou fonction de masse), on note dans la suite \(l(X;\theta) = \ln f(X;\theta)\) et \(l'\), \(l''\) et \(l'''\) désignent respectivement les dérivées première, seconde et troisième par-rapport à \(\theta\).

Théorème Soient \(X_1,\ldots,X_n\) des v.a. IID de densité (ou « fonction de masse ») \(f_0 = f(x;\theta_0) \in \mathcal{F}\) où \(\mathcal{F} = \{f(x;\theta) : \theta \in \Theta\}\) est notre ensemble de densités « candidates ». Si

• \(\hat{\theta}_n \stackrel{P}{\rightarrow} \theta_0\) ;
• \(\hat{\theta}_n\) est tel que \(\sum_i l'(X_i,\hat{\theta}_n) = 0\) (solution des équations de vraisemblance).

Si

1. l'espace des paramètres \(\Theta\) est un sous-ensemble ouvert de \(\mathbb{R}\) ;
2. l'ensemble \(A = \{x : f(x;\theta) > 0\}\) ne dépend pas de \(\theta\) ;
3. \(f(x;\theta)\) est trois fois continument différentiable par-rapport à \(\theta\) pour tout \(x \in A\) ;
4. \(\mathrm{E}_{\theta}\left[l'(X_i;\theta)\right] = 0\) \(\forall\, \theta\) et \(\mathrm{Var}_{\theta}\left[l'(X_i;\theta)\right] = I(\theta)\) où \(0 < I(\theta) < \infty\), \(\forall\, \theta\) ;
5. \(\mathrm{E}_{\theta}\left[l''(X_i;\theta)\right] = -J(\theta)\) où \(0 < J(\theta) < \infty\), \(\forall\, \theta\) ;
6. \(\forall\, \theta\) et \(\delta > 0\), \(|l'''(X_i;t)| \le M(X_i)\) pour \(|t - \theta| \le \delta\) où \(\mathrm{E}_{\theta}\left[M(X_i)\right]< \infty\).

Alors \[\sqrt{n}\left(\hat{\theta}_n - \theta_0\right) \stackrel{L}{\rightarrow} \mathcal{N}\left(0,I(\theta_0)/J^2(\theta_0) \right)\, .\] Quand \(I(\theta_0) = J(\theta_0)\) nous avons : \(\sqrt{n}\left(\hat{\theta}_n - \theta_0\right) \stackrel{L}{\rightarrow} \mathcal{N}\left(0,1/I(\theta_0) \right)\).

Preuve Je ne la reproduit pas ici ; celle donnée en cours suivait le bouquin de Knight.

Supposons que \(X_1,X_2,\ldots,X_n\) sont IID de loi de Cauchy dont nous écrirons la densité : \[f(x;\theta) = \frac{1}{\pi} \frac{1}{1+(x-\theta)^2}\; .\] La fonction de log-vraisemblance est alors : \[\mathcal{l}(\theta) = - n \log \pi - \sum_{i=1}^n \log \left(1+(\theta-X_i)^2\right) \, .\] L'équation de vraisemblance devient : \[\frac{d \mathcal{l}}{d \theta} (\hat{\theta}_n) = \sum_{i=1}^n \frac{2 \, (X_i - \hat{\theta}_n)}{1+(\hat{\theta}_n-X_i)^2} = 0 \, .\]

La loi de Cauchy est l'exemple type d'une loi « pathologique » dans la mesure où aucun de ses moments n'est défini. C'est en fait aussi un exemple assez pathologique pour la méthode du maximum de vraisemblance comme l'a montré Barnett (1966)³. Ici nous allons essentiellement faire ce que les gens ne font presque jamais (et par là comprendre pourquoi les gens ne le font presque jamais), c'est-à-dire vérifier que les hypothèses nécessaires à l'application de nos théorèmes sont remplies. Nous allons ici utiliser le dernier théorème et vérifier les trois « conditions de régularité », (R1), (R2) et (R3). La première se vérifie facilement en calculant les dérivées par-rapport à \(\theta \in \Theta \equiv \mathbb{R}\) et en vérifiant quelles sont finies pour tout \(x\). Avec SymPy (que je charge de façon à ce qu'il interragisse bien avec pylab, c'est-à-dire en gardant un espace de noms séparé pour SymPy) cela donne :

import sympy as sy
x, theta = sy.symbols('x,theta', real=True)
f = 1/(1+(theta-x)**2)/sy.pi

La dérivée première du log de la densité est :

sy.simplify(sy.diff(sy.log(f),theta))

2*(-theta + x)/((theta - x)**2 + 1)

La seconde :

sy.simplify(sy.diff(sy.log(f),theta,2))

2*((theta - x)**2 - 1)/((theta - x)**2 + 1)**2

La troisième :

sy.simplify(sy.diff(sy.log(f),theta,3))

-4*(theta - x)*((theta - x)**2 - 3)/((theta - x)**2 + 1)**3

Où on voit quelles sont toutes définies pour \(x, \theta \in \mathbb{R}\) ce qui vérifie (R1).

Pour vérifier (R2), il faut travailler un peu plus… Tout d'abord nous allons chercher une fonction \(g\) qui domine la norme de la dérivée de la densité (dans un voisinage d'une racine de l'équation de vraisemblance) et qui est intégrable. La dérivé première de la densité est :

sy.simplify(sy.diff(f,theta))

2*(-theta + x)/(pi*((theta - x)**2 + 1)**2)

Si \(\theta_0\) est une solution de l'équation de vraisemblance et si nous considrons un voisinage du type \(V(\theta_0) = (\theta_0-\delta,\theta_0+\delta)\) où \(\delta > 0\), alors, pour \(|x-\theta_0| < \delta\), \(|-\theta +x| \le 2 \delta\). Pour \(x-\theta_0 \ge \delta\) nous voyons que \(\partial f(x;\theta) / \partial \theta \le 2 (x + \delta - \theta_0) / (\pi ((x - \delta - \theta_0)^2+1)^2)\) ; ainsi : \[g(x) = \left\{\begin{array}{l l} \frac{2 (-x + \delta + \theta_0)}{\pi ((x - \delta - \theta_0)^2+1)^2} & \quad \text{si} \quad -x+\theta_0 \le -\delta \\ 2 \delta & \quad \text{si} \quad |x-\theta_0| < \delta \\ \frac{2 (x + \delta - \theta_0)}{\pi ((x - \delta - \theta_0)^2+1)^2} & \quad \text{si} \quad x-\theta_0 \ge \delta \end{array} \right.\] domine la valeur absolue de la dérivée première pour tout \(\theta \in V(\theta_0)\) et est clairement intégrable sur \(\mathbb{R}\).

Quand on fait ce genre de calcul, il peut être prudent de faire quelques graphes afin de vérifier que notre fonction dominante, domine bien ce qu'elle est sensée dominer… On pourra le faire ici de la façon suivante (les graphes ne sont pas montrés dans ce document) :

xx = linspace(-5,5,501)

def abs_df(x,theta):
    """Renvoie la valeur absolue de la derivee premiere de la densite"""
    return abs(2*(-theta + x)/(pi*((theta - x)**2 + 1)**2))

def g(x,delta,theta0=0):
    """Renvoie une fonction dominant abs_df sur un voisinage de longueur delta autour de theta0"""
    if abs(x-theta0) < delta:
        return 2*delta
    else:
        return 2*(abs(x-theta0)+delta)/pi/((abs(x-theta0)-delta)**2+1)**2

plot(xx,[abs_df(x,0) for x in xx],color='black')
plot(xx,[abs_df(x,0.5) for x in xx],color='black')
plot(xx,[abs_df(x,-0.5) for x in xx],color='black')
plot(xx,[g(x,0.5) for x in xx],color='red')

La dérivé seconde de la densité est :

sy.simplify(sy.diff(f,theta,2))

2*(3*(theta - x)**2 - 1)/(pi*((theta - x)**2 + 1)**3)

On voit alors que : \[h(x) = \left\{\begin{array}{l l} \frac{6 (|x-\theta_0| + \delta)^2}{\pi ((|x - \theta_0| - \delta )^2+1)^3} & \quad \text{si} \quad |x-\theta_0| \ge 1/\sqrt{3} \\ 2 / \pi & \quad \text{si} \quad |x-\theta_0| < 1/\sqrt{3} \end{array} \right.\] domine la valeur absolue de la dérivée seconde et est intégrable.

On peut encore vérifier graphiquement qu'on ne sait pas trompé avec :

def abs_ddf(x,theta):
    """Renvoie la valeur absolue de la derivee seconde de la densite"""
    return abs(2*(3*(theta - x)**2 - 1)/(pi*((theta - x)**2 + 1)**3))

def h(x,delta,theta0=0):
    """Renvoie une fonction dominant abs_ddf sur un voisinage de longueur delta autour de theta0"""
    if abs(x-theta0) < 1/sqrt(3):
        return 2/pi
    else:
        return 6*(abs(x-theta0)+delta)**2/pi/((abs(x-theta0)-delta)**2+1)**3

plot(xx,[abs_ddf(x,0) for x in xx],color='black')
plot(xx,[abs_ddf(x,0.5) for x in xx],color='black')
plot(xx,[abs_ddf(x,-0.5) for x in xx],color='black')
plot(xx,[h(x,0.5) for x in xx],color='red')
plot(xx,[h(x,sqrt(3)) for x in xx],color='blue')

La dérivée troisième du log de la densité est :

sy.simplify(sy.diff(sy.log(f),theta,3))

-4*(theta - x)*((theta - x)**2 - 3)/((theta - x)**2 + 1)**3

Si nous choisissons \(V(\theta_0) = (\theta_0 - \sqrt{3},\theta_0 + \sqrt{3})\), choix compatible avec les deux cas précédents, alors : \[H(x) = \left\{\begin{array}{l l} \frac{4 (|x-\theta_0| + \sqrt{3})}{((|x - \theta_0| - \sqrt{3})^2+1)^2} & \quad \text{si} \quad |x-\theta_0| \ge \sqrt{3} \\ 12 \sqrt{3} & \quad \text{si} \quad |x-\theta_0| < \sqrt{3} \end{array} \right.\] domine la valeur absolue de la dérivée troisième du log de la densité. Il nous reste à montrer que l'espérance de \(H(x)\) par-rapport à toute densité \((\pi (1+(\theta - x)^2))^{-1}\) où \(\theta \in (\theta_0 - \sqrt{3},\theta_0 + \sqrt{3})\) existe. Remarquez, pour les quelques calculs qui restent, que nous avons une invariance par translation (c'est toujours la différence \(x - \theta_0\) qui apparaît) et nous pouvons ainsi poser \(\theta_0 = 0\). Alors le morceau central de \(H(x)\) (pour \(x < \sqrt{3}\)) ne pose pas de problème d'intégration par-rapport à la densité. Il nous faut considérer les deux morceaux (principaux) symétriques, celui de droite étant : \[\int_{\sqrt{3}}^{\infty} \frac{4 (x + \sqrt{3})}{((x - \sqrt{3})^2+1)^2} \frac{1}{\pi (1+(\theta - x)^2)} dx \, ,\] où \(\theta \in (-\sqrt{3},\sqrt{3})\). Mais on a clairement : \[\int_{\sqrt{3}}^{\infty} \frac{4 (x + \sqrt{3})}{((x - \sqrt{3})^2+1)^2} \frac{1}{\pi (1+(\theta - x)^2)} dx \le \int_{\sqrt{3}}^{\infty} \frac{4 (x + \sqrt{3})}{((x - \sqrt{3})^2+1)^2} \frac{1}{\pi (1+(\sqrt{3} - x)^2)} dx \le \frac{4}{\pi} \int_{\sqrt{3}}^{\infty} \frac{x + \sqrt{3}}{(x - \sqrt{3})^6+1} dx\] et cette dernière intégrale vaut :

x = sy.symbols('x',real=True)
sy.simplify(4/sy.pi*sy.integrate((x+sy.sqrt(3))/(1+(x-sy.sqrt(3))**6),(x,sy.sqrt(3),sy.oo)))

28*sqrt(3)/9

Et on vérifie graphiquement que cela colle avec :

def abs_dddlogf(x,theta):
    """Renvoie la valeur absolue de la derivee troisieme de la log densite"""
    return abs(-4*(theta - x)*((theta - x)**2 - 3)/((theta - x)**2 + 1)**3)

def H(x,theta0=0):
    """Renvoie une fonction dominant abs_dddlogf sur un voisinage de longueur sqrt(3) autour de theta0"""
    if abs(x-theta0) < sqrt(3):
        return 12*sqrt(3)
    else:
        return 4*(abs(x-theta0)+sqrt(3))/((abs(x-theta0)-sqrt(3))**2+1)**2

plot(xx,[abs_dddlogf(x,0) for x in xx],color='black')
plot(xx,[abs_dddlogf(x,1) for x in xx],color='black')
plot(xx,[abs_dddlogf(x,sqrt(3)) for x in xx],color='black')
plot(xx,[abs_dddlogf(x,-1) for x in xx],color='black')
plot(xx,[abs_dddlogf(x,-sqrt(3)) for x in xx],color='black')
plot(xx,[H(x) for x in xx],color='red')

Ce qui achève la vérification de (R2).

Grâce à SymPy, vérifier (R3) est un jeu d'enfant, puisque nous obtenons l'expression analytique de l'espérance du carré de la dérivé première (par-rapport au paramètre) du log de la densité avec :

x = sy.symbols('x',real=True)
sy.simplify(sy.integrate(sy.diff(sy.log(f),theta)**2*f,(x,-sy.oo,sy.oo)))

1/2

De façon un peu formelle, suivant le bouquin de Lejeune (p 136-137), on définit :

Définition Soit \(X_1,X_2,\ldots,X_n\) un échantillon issu d'une densité (ou fonction de masse) \(f(x;\theta_0)\) où \(\theta_0 \in \Theta\) est un paramètre inconnu de dimension 1. On appelle procédure d'intervalle de confiance de niveau \(\alpha\) tout couple de statistiques⁵ \((T_1,T_2)\) tel que, quelque soit \(\theta_0 \in \Theta\), on ait : \[\mathrm{Prob}_{\theta_0}\left(T_1 \le \theta_0 \le T_2\right) \ge \alpha \, .\]

Définition Dans le contexte de la définition précédente, soit \(x_1,x_2,\ldots,x_n\) une réalisation de \(X_1,X_2,\ldots,X_n\) conduisant à la réalisation \((t_1,t_2)\) de \((T_1,T_2)\). Alors l'intervalle \([t_1,t_2]\) est appelé intervalle de confiance de niveau \(\alpha\) pour \(\theta_0\) et l'on note : \[IC_{\alpha}(\theta_0) = [t_1,t_2] \, .\]

Remarque L'intervalle de confiance n'est rien d'autre que l'application numérique de la procédure suite à la réalisation de l'échantillon.

La densité exponentielle « complète » est \(p(t) = \nu \exp ( - t \nu)\), la fonction de répartition est donc \[ P(t) = 1 - \exp ( - t \nu) \] et \[ \textrm{Pr}\{T > \delta\} = \exp ( - \delta \nu) \; . \] La densité de probabilité supposée pour le problème est donc : \[ p(t \mid t \ge \delta) = \nu \exp \left( - (t-\delta) \nu\right)\; . \] La log vraisemblance du problème est ainsi : \[ l_n(\nu) = N \, \log \nu - \nu \sum_i (t_i-\delta) \; , \] où \(N=20\) et \(i=1,\ldots,20\). La fonction score (dérivée première) est : \[ l_n'(\nu) = N/\nu - \sum_i (t_i - \delta) \] et la dérivée seconde est \[ l_n''(\nu) = -N/\nu^2\; . \] Cette dernière est toujours négative (tant que \(N>0\), ce qui est le cas ici) et \(l_n'\) est donc concave et admet un maximum unique donné par : \[ \hat{\nu} = N / \sum_i (t_i - \delta) \; . \] Numériquement il vient :

nu_chapeau = 1/(mean(intervalles)-delta)
nu_chapeau

0.77315602288541818

Le développement précédent nous donne : \[ \widehat{es}_{\hat{\nu}} = \hat{\nu} / \sqrt{N} \; , \] soit :

es_nu_chapeau = nu_chapeau/sqrt(n_intervalles)
es_nu_chapeau

0.17288294243851782

D'où l'intervalle de confiance à 95% avec la méthode de Wald :

np.round(nu_chapeau + 1.96*array([-1,1])*es_nu_chapeau,decimals=2)

array([ 0.43,  1.11])

Nous traçons le graphe de , \(l_n(\nu)\), pour \(\nu \in [-3*\sigma_{\hat{\nu}} + \hat{\nu},-3*\sigma_{\hat{\nu}} + \hat{\nu}]\). Nous y rajoutons celui de l'approximation quadratique que nous venons de faire (implicitement) pour obtenir votre intervalle de confiance : \(y = l_n(\hat{\nu}) - (\nu - \hat{\nu})^2/(2\, \widehat{es}_n^2)\).

def log_vrai(nu): return n_intervalles*((delta-mean(intervalles))*nu+log(nu))
nn = linspace(-3,3,101)*es_nu_chapeau + nu_chapeau
plot(nn,[log_vrai(nu) for nu in nn],lw=2,color='black')
xlabel('$\\nu$')
ylabel('Log vraisemblance')
y_max = log_vrai(nu_chapeau)
y_min = min(log_vrai(nn[0]),log_vrai(nn[-1]))
plot(nn,[y_max-(nu-nu_chapeau)**2/2/es_nu_chapeau**2 for nu in nn],lw=2,color='red')
vlines(-1.96*es_nu_chapeau + nu_chapeau,y_min,y_max,color='black',linestyle='dashed')
vlines(1.96*es_nu_chapeau + nu_chapeau,y_min,y_max,color='black',linestyle='dashed')

<matplotlib.collections.LineCollection object at 0x4b615d0>

savefig('img/coursIntroStats2014/comparaison-logvrai-et-quadratique-compteur-temps-mort.png')
close()
'img/coursIntroStats2014/comparaison-logvrai-et-quadratique-compteur-temps-mort.png'

On voit que l'approximation quadratique commence à être « assez loin » de la log vraisemblance même dans le domaine défini par l'intervalle de confiance. Qualitativement, la convergence en loi de l'EMV vers une normale peut s'interpréter en disant que, si la taille de l'échantillon est suffisamment grande, le graphe de la log vraisemblance doit être très proche d'une parabole retournée. Le fait que le graphe ci-dessus s'éloigne « assez vite » de la parabole qui l'approche en son maximum, suggère que la taille de l'échantillon est trop petite pour que l'intervalle obtenu avec la méthode de Wald soit fiable.

Numériquement, \(k_{0,05}\) vaut :

from scipy.stats import chi2
k_0p05 = chi2.ppf(1-0.05,1)
k_0p05

3.8414588206941236

Nous allons donc chercher les abscisses des points d'intersection entre \(l_n(\nu)\) et une droite horizontale d'ordonnée \(l_n(\hat{\nu}) - k_{0,05}/2\). Pour trouver ces abscisses, nous devons trouver les racines de la fonction (de \(\nu\)) : \[l_n(\hat{\nu}) - l_n(\nu) - k_{0,05}/2 \] pour \(\nu \in [0,2;\hat{\nu}]\) d'une part et pour \(\nu \in [\hat{\nu};1,2]\) d'autre part. Nous allons faire cela avec la méthode de Brent mise en œuvre par la fonction brentq de scipy :

from scipy.optimize import brentq
def cible(nu): return log_vrai(nu_chapeau) - log_vrai(nu) - k_0p05/2.
IC_lr_g = brentq(cible,0.2,nu_chapeau)
IC_lr_d = brentq(cible,nu_chapeau,1.2)
[IC_lr_g,IC_lr_d]
 

[0.4818964626882697, 1.163206521355929]

plot(nn,[log_vrai(nu) for nu in nn],lw=2,color='black')
xlabel('$\\nu$')
ylabel('Log vraisemblance')
hlines(log_vrai(nu_chapeau)-k_0p05/2.,0.2,1.4,color='red',lw=2)
vlines(-1.96*es_nu_chapeau + nu_chapeau,y_min,y_max,color='black',linestyle='dashed')
vlines(1.96*es_nu_chapeau + nu_chapeau,y_min,y_max,color='black',linestyle='dashed')
vlines(IC_lr_g,y_min,y_max,color='red',linestyle='dashed')
vlines(IC_lr_d,y_min,y_max,color='red',linestyle='dashed')

<matplotlib.collections.LineCollection object at 0x4c61ed0>

savefig('img/coursIntroStats2014/logvrai-et-logratio-compteur-temps-mort.png')
close()
'img/coursIntroStats2014/logvrai-et-logratio-compteur-temps-mort.png'

On voit que l'intervalle de la méthode du rapport de vraisemblance est « décalé » vers la droite par-rapport à l'intervalle obtenu avec l'approximation normale. Ceci est cohérent avec la différence entre l'approximation quadratique et la fonction de log vraisemblance.

Vous prendrez la valeur de \(\hat{\nu}\) obtenue comme la vraie valeur. Vous simulerez de 500 à 1000 échantillons de taille 20 avec cette valeur. Pour chaque échantillon simulé, vous calculerez les \(IC_{0,95}\) avec la méthode de Wald et avec le méthode du rapport de vraisemblance. Vous obtiendrez ainsi les probabilités de recouvrement empiriques de vos intervalles. Conclusion ?

Pour ce type d'exercice ou d'étude basés sur la répétition d'une même analyse sur une collection de réplique, je vous recommande de commencer par définir de courtes fonctions qui effectuent une partie bien précise de l'analyse demandée. Ici, pour chaque réplique (échantillon de 20 observations simulées) nous allons construire un \(IC_{0,95}\) avec la méthode de Wald et un autre avec la méthode du rapport de vraisemblance. Il me semble donc raisonnable de commencer par définir deux fonctions, une qui, pour une réplique donnée, renvoie l'\(IC_{0,95}\) avec la méthode de Wald, l'autre qui renvoie l'\(IC_{0,95}\) obtenu avec la méthode du rapport de vraisemblance. La première des deux fonctions est simple à définir, il suffit de retrouver les quelques lignes de code utilisées ci-dessus pour calculer l'intervalle sur les données originales et de remplacer dans le code les données originales par l'argument de la fonction :

def IC_wald(data):
    nu_chapeau = 1/(mean(data)-delta)
    es_nu_chapeau = nu_chapeau/sqrt(n_intervalles)
    return [nu_chapeau - 1.96*es_nu_chapeau,nu_chapeau + 1.96*es_nu_chapeau]

Remarquez que la fonction retourne une liste dont les deux éléments sont les bornes de notre intervalle de confiance. Quand nous avions fait les calculs pas à pas sur les données originales, nous avions les bornes contenues dans un vecteur numpy. La liste, plus classique, choisie ici rend le traitement ultérieur un peu plus simple, c'est pourquoi elle a été choisie.

Quand on développe de courtes fonctions comme IC_wald, on a tout intérêt à effectuer de petits tests pour s'assurer qu'aucune erreur margeur n'a été commise. Ici, la première à faire est de vérifier que notre fonction nous redonne l'intervalle calculé précédemment « à la main » quand on l'applique aux mêmes données, les données originales :

IC_wald(intervalles)

[0.43430545570592327, 1.112006590064913]

Nous continuons en définissant une fonction qui calcule un \(IC_{0,95}\) par la méthode du rapport de vraisemblance. Ici, la tâche est un peu plus délicate car le calcul nécessite deux appels à brentq, appels pour lesquels une fonction cible doit être fournie. Cette fonction dépend des données et va donc devoir être redéfinie pour chaque nouveau jeu de données ; nous allons donc devoir la définir comme une fonction locale. Enfin, les bornes entre lesquelles la racine est cherchée doivent être fournies lors des deux appels à brentq. La borne droite, lors du premier appel (qui est la borne gauche lors du second est connue) est connue lors de l'appel puisqu'elle correspond à l'EMV. Le problème est le choix de la borne gauche (et celui de la borne droite lors du second appel). Nous avons précédemment résolu ce problème en examinant le graphe de la fonction, mais nous ne voulons pas utiliser cette approche ici puisque nous voulons une génération automatique de l'\(IC_{0,95}\). Une première stratégie consiste à prendre des bornes assez éloignées de l'EMV. En considérant que l'approximation quadratique n'est pas trop mauvaise, on se dit que fixer les bornes à un plus que 2 fois l'erreur type (\(\widehat{es}\)) ne devrait pas être trop mal. Dans un premier temps, j'avais décidé d'être « généreux » et j'avais fixer les bornes à \(5 \times \widehat{es}\) de part et d'autre de l'EMV. Là, les premiers essais m'ont générés des erreurs car le log d'un nombre négatif était demandé : en clair, l'EMV était « trop » proche de 0 et \(\widehat{es}\) était « trop » large ce qui donnait \(\hat{\nu} - 5 \times \widehat{es}\). J'ai donc décider de fixer la limite gauche lors du premier appel à brentq au maximum de \(\hat{\nu} - 5 \times \widehat{es}\) et d'un nombre suffisamment petit que j'ai fixé à 0,1. D'où la fonction :

def IC_rv(data):
    nu_chapeau = 1/(mean(data)-delta)
    def log_vrai(nu): return n_intervalles*((delta-mean(data))*nu+log(nu))
    def cible(nu): return log_vrai(nu_chapeau) - log_vrai(nu) - k_0p05/2.
    es_nu_chapeau = nu_chapeau/sqrt(n_intervalles)
    IC_lr_g = brentq(cible,max(nu_chapeau-5*es_nu_chapeau,0.1),nu_chapeau)
    IC_lr_d = brentq(cible,nu_chapeau,nu_chapeau+5*es_nu_chapeau)
    return [IC_lr_g,IC_lr_d]

Là encore, nous vérifions qu'appliquée aux données originales, notre fonction nous redonne l'\(IC_{0,95}\) obtenue précédemment « à la main » :

IC_rv(intervalles)

[0.4818964626882665, 1.1632065213559313]

Nous pouvons alors effectuer notre série de simulations. Mais avant de nous lancer dans une simulation à taille réelle, il vaut mieux en faire une avec peu de répliques (disons 10) et vérifier que rien de franchement pathologique n'apparaît. Après seulement, on peut se lancer avec beaucoup de répliques. Ici j'ai en fait constaté que les calculs aller tellement vite que j'ai fini par faire 10000 répliques au lieu de 1000 :

wald = []
rv = []
nu_sim = 1/(mean(intervalles)-delta)
for i in range(10000):
    echantillon = expon.ppf(uniform.rvs(loc=expon.cdf(delta,scale=1/nu_sim),scale=1-expon.cdf(delta,scale=1/nu_sim),size=n_intervalles),scale=1/nu_sim)
    echantillon = np.round(echantillon,decimals=3)
    wald.append(IC_wald(echantillon))
    rv.append(IC_rv(echantillon))

Nous voulons maintenant, pour chacun des intervalles obtenus, tester si la vraie valeur (nu_sim) et, oui ou non, contenue dedans, puis compter combien de fois, elle l'est. Pour les intervalles calculés avec la méthode de Wald nous obtenons :

sum([ic[0] <= nu_sim <= ic[1] for ic in wald])

9538

Pour ceux calculés avec la méthode du rapport de vraisemblance nous avons :

sum([ic[0] <= nu_sim <= ic[1] for ic in rv])

9510

Dans ce cas, l'espérance du nombre de oui est 10000 x 0,95 = 9500 ; nous faisons donc un peu mieux avec la seconde méthode, mais la première ne donne pas des résultats catastrophiques.

Cet exemple est clairement « un peu » artificiel (du point de vue méthodologique). En effet, comme l'EMV est une fonction de la moyenne empirique, nous aurions mieux fait de procéder comme dans la section 4.4.2. Nous aurions même pu décider de reparamétriser la loi des données en travaillant avec le temps caractéristique de l'exponentielle, plutôt que sa fréquence caractéristique.

• Question 2.1 : vous écrirez la log-vraisemblance, la fonction score (dérivée de la log-vraisemblance par-rapport au paramètre) et l'information observée (opposée de la dérivée seconde de la log-vraisemblance par-rapport au paramètre) pour un échantillon de taille N, IID, suivant une loi de Rayleigh. Vous en déduirez l'expression analytique de l'estimateur du maximum de vraisemblance ;
• Question 2.2 : pour chacun des trois jeux de données vous trouverez l'estimateur du maximum de vraisemblance du paramètre θ de la loi de Rayleigh ainsi qu'un intervalle de confiance à 95% avec la méthode de Wald ;
• Question 2.3 : pour chacun des trois jeux de données vous construirez le graphe de la fonction de log-vraisemblance au voisinage de son maximum (sur un domaine de longueur comparable à celle de l'intervalle de confiance construit, par exemple, 3 fois l'erreur type) auquel vous superposerez l'approximation quadratique : \[y = l(\hat{θ}) - (\hat{θ} - θ)^2 / (2 \widehat{es}^2)\; ; \]
• Question 2.4 : vous obtiendrez un intervalle de confiance à 95% avec la méthode du rapport de vraisemblance pour chacun des trois jeux de données ;
• Question 2.5 : pour chacun des trois jeux de données vous construirez un histogramme des observations (en choisissant la largeur des classes par validation croisée) auquel vous superposerez la loi de Rayleigh ajustée. Le modèle ajusté vous semble-t-il satisfaisant ?
• Question 2.6 : pour chacun des trois jeux, vous construirez le graphe de la fonction de répartition empirique, dont le le code se trouve au début du cours, à laquelle vous rajouterez une bande de confiance à 95% (le code se trouve peu après celui définissant la construction du graphe de la fonction de répartition empirique) avant de rajouter la fonction de répartition ajustée ;
• Question 2.7 : vous construirez un graphe représentant votre paramètre estimé \(\hat{θ}\) en fonction de la racine carrée de la distance linéaire séparant les marqueurs. Qu'en concluez-vous sachant que l'écart type de la distance euclidienne entre l'origine et la fin d'une marche aléatoire croît comme la racine du nombre de pas quelle comporte ?
• Question 2.8 : pour un des jeux de données, vous étudierez la probabilité de recouvrement d'un intervalle de confiance à 95% obtenu avec la méthode de Wald par bootstrap paramétrique (comme nous l'avons fait dans l'exercice de la section 4.7.4).

Le module scipy.stats fournit une collection de méthodes qui vous permettent de calculer la densité, la fonction de répartition et de générer des nombres (pseudo-)aléatoires suivant une loi de Rayleigh.