Examen du 18 juin 2015 : une question de temps au kilomètre pour les baleines

Les données

Vous allez étudier un jeu de données relatif à des baleines boréales, Balaena mysticetus. Ici des repérages d'individus identifiés ont été effectués successivement à deux positions et la vitesse de déplacement (en km/h) en a été déduite. Le jeux de données contient l'inverse de cette vitesse, soit le temps (en heures) nécessaire pour parcourir 1 km. Après avoir débuté votre session Python suivant la « procédure habituelle » :

In [1]:
%pylab
from __future__ import print_function, division, unicode_literals, absolute_import
%matplotlib inline

Using matplotlib backend: Qt4Agg
Populating the interactive namespace from numpy and matplotlib

Vous chargerez les données de la façon suivante :

In [2]:
from urllib2 import urlopen
u = urlopen("http://xtof.disque.math.cnrs.fr/cours/whales.txt")
baleines = [float(x) for x in u.readlines()]
u.close()

Format des réponses aux questions suivantes

Vous remplirez un notebook IPython avec vos lignes de codes et, si possible, quelques lignes de textes de temps en temps pour expliquer ce que vous faites. Vous sauverez votre notebook que vous pourrez intituler de vos nom et prénom et, si possible, vous l'exporterez au format html. À la fin de l'examen, vous m'enverrez au moins votre notebook, c'est-à-dire un fichier au format .ipynb (soit un fichier du genre nom_prenom.ipynb) et, si possible, la même chose au format html. Vous enverrez le ou les fichiers à l'adresse : christophe.pouzat@parisdescartes.fr. En cas de problème, j'aurais aussi une clé USB pour récupérer votre travail.

Histogramme (Question 1 -- poids 1)

Vous commencerez par construire un histogramme des données. Vous essaierez plusieurs nombres de classes / bins et en choisirez une qui vous semble « raisonnable » -- je ne vous demande pas ici d'utiliser les méthodes sophistiqués présentées lors du cours de l'an dernier, juste un choix indicatif « à l'œil ». Nous avons déjà construit un histogramme avec la fonction hist lors du cours du 7 mai dernier. N'oubliez pas d'indiquer ce qui est montré sur les axes (par exemple Temps (heure)).

Réponse

À ce stade, les données baleines constituent une liste. Pour faciliter les manipulations qui vont suivre, je vais les convertir en vecteur (array) numpy :

In [3]:
baleines = array(baleines)

Je construit alors l'histogramme que je normalise en spécifiant normed = True (ici, vous êtes libre de normaliser ou non, cela va juste jouer sur la labélisation de votre ordonnée).

In [9]:
baleines_hist = hist(baleines,bins=50,normed=True,color='black')
xlabel('Temps (heure)')
ylabel(u'Densité estimée')
Out[9]:
<matplotlib.text.Text at 0x7fd246245b90>

J'ai construit ici des histogrammes avec 10, 25, 50, 100 et 200 bins / classes. Pour moi, 10, c'est trop peu et 200, c'est trop, entre les deux, c'est OK.

Loi Gamma et méthode des moments (Question 2 -- poids 2)

L'article d'où je tire ces données, présente une analyse basée sur une loi gamma. Vous pourrez trouver de nombreuses informations importantes sur cette loi sur la page qui lui est consacrée sur wikipedia (en version anglaise sous le titre Gamma distribution, la française étant beaucoup plus pauvre). La loi gamma est disponible dans le sous-module scipy.stats. Python emploie deux paramètres pour cette loi:

  • a, correspondant au paramètre k de la page wikipedia ;
  • scale, correspondant au paramètre $\theta$.

C'est-à-dire que la densité de probabilité que wikipedia écrit : $x^{k-1}e^{-x/\theta}/\Gamma(k)/\theta^k$ ; la documentation Python l'écrit : $x^{a-1}e^{-x/\mathrm{scale}}/\Gamma(a)/\mathrm{scale}^a$.

Méthode des moments

La méthode des moments consiste en (dans le cas qui nous occupe) :

  1. Trouver des expressions de l'espérance et de la variance d'une variable aléatoire de loi gamma de paramètres a et scale ($X \sim \mathcal{G}(a, \mathrm{scale})$) en fonction de ces paramètres : $\mathrm{E}X = f(a,\mathrm{scale})$ et $\mathrm{Var}X = g(a,\mathrm{scale})$ ;
  2. Obtenir des estimations empiriques : $\overline{X} \approx \mathrm{E}X$ et $S^2 = \sum_{i=1}^n (X_i-\overline{X})^2/(n-1) \approx \mathrm{Var}X$ ;
  3. Trouver $\tilde{a}$ et $\widetilde{\mathrm{scale}}$, les estimateurs de la méthode des moments, tels que $\overline{X} = f(\tilde{a},\widetilde{\mathrm{scale}})$ et $S^2 = g(\tilde{a},\widetilde{\mathrm{scale}})$.

Votre tâche

  1. Vous trouverez $\tilde{a}$ et $\widetilde{\mathrm{scale}}$ pour les données ci-dessus (la page wikipedia vous fournira les expression $f$ et $g$).
  2. Vous définirez une fonction python qui prendra comme argument un jeu de données et retourne un tuple (une paire) constitué de $\tilde{a}$ et $\widetilde{\mathrm{scale}}$.

Réponses

Sur la page wikipedia nous voyons que : $\mathrm{E}X = a \times \mathrm{scale}$ et que $\mathrm{Var}X = a \times \mathrm{scale}^2$. Nous devons donc résoudre : $$\overline{X} = a \times \mathrm{scale}$$ et $$S^2 = a \times \mathrm{scale}^2$$ pour les deux paramètres, ce qui nous donne : $$\widetilde{\mathrm{scale}} = S^2/\overline{X}$$ et $$\tilde{a} = \overline{X}^2/S^2 \, .$$ Soit en Python :

In [10]:
baleines_moy = mean(baleines)
baleines_var = var(baleines)
a_tilde,scale_tilde = (baleines_moy**2/baleines_var,baleines_var/baleines_moy)
(a_tilde,scale_tilde)
Out[10]:
(0.79917407513851624, 0.7582816641218163)

Nous définissons maintenant une courte fonction faisant tout cela :

In [11]:
def methode_moments(donnees):
    x_bar = mean(donnees)
    S2 = var(donnees)
    a = S2/x_bar
    scale = x_bar**2/S2
    return (a, scale)

Nous la testons avec le cas que nous avons juste traité :

In [12]:
methode_moments(baleines)
Out[12]:
(0.7582816641218163, 0.79917407513851624)

Estimation des paramètres par la méthode du maximum de vraisemblance (Question 3 -- poids 3)

Vous lirez la section Maximum likelihood estimation de la page wikipedia (sur la Gamma distribution). Vous coderez une ou des fonctions Python vous permettant d'obtenir l'EMV (estimateur du maximum de vraisemblance) en utilisant la méthode de Newton pour le paramètre a. Vous aurez pour cela besoin des fonctions digamma et polygamma du sous-module scipy.special. Nous avons fait quelque chose de très semblable (méthode de Newton, etc) dans le cours du 9 avril dernier. Vous ferez attention à vérifier que votre algorithme de Newton s'arrête bien parce-qu'une solution a été trouvée et pas par-ce que le nombre maximal d'itérations a été atteint.

De même que pour la méthode des moments, vous définirez une fonction Python prenant en entrée, un jeu de données et retournant en sortie un triplé dont les éléments seront : $\hat{a}$, $\widehat{\mathrm{scale}}$ et le nombre d'itérations de la méthode de Newton. Ici $\hat{a}$ et $\widehat{\mathrm{scale}}$ désignent les estimateur obtenus par la méthode du maximum de vraisemblance.

Réponse

La page wikipedia nous donne : $$\widehat{\mathrm{scale}}(a) = \frac{1}{a n} \sum_{i=1}^n x_i \, ,$$ c'est-à-dire $$\widehat{\mathrm{scale}}(a) = \overline{X}/a \, .$$ Le pas de la méthode de Newton fait intervenir une quantité $s = \log \left(\sum_{i=1}^n x_i / n\right) - 1/n \sum_{i=1}^n \log x_i$, c'est-à-dire, la différence entre le log de la moyenne et la moyenne des log. Commençons par calculer s pour notre échantillon :

In [13]:
s = log(mean(baleines)) - mean(log(baleines))
s
Out[13]:
0.34503126893469949

Le numérateur du second terme du membre de droite du pas de Newton fait intervenir la fonction digamma alors que le dénominateur du même terme fait intervenir sa dérivée -- que nous pouvons calculer grâce à la fonction polygamma. Suivant le cours du 9 avril, nous définissons une fonction effectuant une séquence de pas de Newton jusqu'à ce qu'une valeur suffisamment proche de la racine soit trouvée ou jusqu'à ce que 100 itérations aient été effectuées (c'est juste un copié-collé de la fonstion du cours) :

In [14]:
def newton(devinette_initiale,
           cible,
           cible_derivee,
           tolerance=1e-6,
           iter_max=100):
    pos = devinette_initiale
    valeur = cible(pos)
    idx = 0
    while idx <= iter_max and abs(valeur) > tolerance :
        pos -= valeur/cible_derivee(pos)
        idx += 1
        valeur = cible(pos)
    return (pos,valeur,idx)

Il nous reste à définir la fonction cible, second argument de la fonction newton que nous venons de définir (mais aussi numérateur du second terme du membre de droite de l'équation définissant un pas de newton sur la page wikipedia) et la fonction cible_derive (dénominateur du second terme). Nous effectuons cela simplement, ou presque avec :

In [15]:
from scipy.special import digamma, polygamma
def num_fct(a,s=s):
    return log(a) - digamma(a) - s

def denom_fct(a):
    return 1/a - polygamma(1,a)

Il nous faut encore trouver la devinette initiale pour a, mais il suffit encore de suivre la page wikipedia pour avoir :

In [20]:
a_0 = (3-s+sqrt((s-3)**2+24*s))/12/s
a_0
Out[20]:
1.5868786962409753

Nous obtenons à présent nos estimateurs :

In [19]:
a_hat, cible_val, nb_iter = newton(a_0,num_fct,denom_fct)
(a_hat, cible_val, nb_iter)
Out[19]:
(1.5954085788795105, 3.7776948236256658e-10, 2)

Nous voyons que deux itérations ont été suffisantes pour converger. Nous obenons maintenant scale_hat :

In [21]:
scale_hat = baleines_moy/a_hat
scale_hat
Out[21]:
0.37983940643258518

Nous pouvons maintenant définir une fonction qui fait tout ce travail d'un coup :

In [22]:
def methode_mv(donnees):
    x_bar = mean(donnees)
    s = log(x_bar) - mean(log(donnees))
    from scipy.special import digamma, polygamma
    def num_fct(a):
        return log(a) - digamma(a) - s
    def denom_fct(a):
        return 1/a - polygamma(1,a)
    a_0 = (3-s+sqrt((s-3)**2+24*s))/12/s
    a_hat, cible_val, nb_iter = newton(a_0,num_fct,denom_fct)
    scale_hat = x_bar/a_hat
    return (a_hat,scale_hat, nb_iter)

Nous vérifions que nous obtenons avec cette fonction les résultats que nous avons déjà obtenus pas à pas :

In [23]:
methode_mv(baleines)
Out[23]:
(1.5954085788795105, 0.37983940643258518, 2)
In [24]:
(a_hat,scale_hat, nb_iter)
Out[24]:
(1.5954085788795105, 0.37983940643258518, 2)

Histogramme et loi ajustée (Question 4 -- poids 2)

Vous construirez une figure avec l'histogramme précédent et les deux densités obtenues avec les deux méthodes d'estimations. Commentez.

In [26]:
from scipy.stats import gamma
tt = linspace(0,4.5,201)
yy_tilde = gamma.pdf(tt,a=a_tilde, scale=scale_tilde)
yy_hat = gamma.pdf(tt,a=a_hat, scale=scale_hat)
baleines_hist = hist(baleines,bins=50,normed=True,color='black')
xlabel('Temps (heure)')
ylabel(u'Densité estimée')
plt.plot(tt,yy_tilde,color='blue',lw=2)
plt.plot(tt,yy_hat,color='red',lw=2)
Out[26]:
[<matplotlib.lines.Line2D at 0x7fd23db84410>]

La densité obtenue par l'EMV est « moins mauvaise » -- dans la mesure où elle est plus grande là ou l'hitogramme a de grandes valeurs -- mais aucune des deux n'a l'air vraiment satisfaisante.

Pour quantifier la différence de qualité des deux estimations, nous pouvons calculer leurs log-vraisemblances respectives :

In [28]:
(sum(log(gamma.pdf(baleines,a=a_tilde,scale=scale_tilde))),
 sum(log(gamma.pdf(baleines,a=a_hat,scale=scale_hat))))
Out[28]:
(-117.81409799069898, -92.723789284357466)

La log vraisemblance obtenue avec l'EMV est bien plus grande.

Estimation des loi d'échantillonage par méthode du bootstrap (Question 5 -- poids 4)

Dans le cours nous avons utilisé plusieurs fois un bootstrap paramétrique ; nous allons ici faire quelque chose d'encore plus simple, un bootstrap de base. L'idée est la suivante :

  • nous disposons d'un échantillon $x_1,\ldots,x_n$ d'observations de variables aléatoires $X_1,\ldots,X_n$ supposées IID et de loi inconnue ;
  • nous avons des statistiques, c'est-à-dire des fonctions définies sur cet échantillon, $S(x_1,\ldots,x_n)$ -- ici $\tilde{a}$, $\widetilde{\mathrm{scale}}$, $\hat{a}$ et $\widehat{\mathrm{scale}}$ sont des exemples de statistiques -- ;
  • nous souhaitons estimer les lois des $S(X_1,\ldots,X_n)$ qui, elles, sont des variables aléatoires ;
  • le bootstrap obtient une telle estimation en ré-échantillonnant les données.

Le ré-échantillonage des données, $x_1,\ldots,x_n$, consiste en la génération d'un « nouvel » échantillon, $y_1,\ldots,y_n$, (dit de bootstrap) où les $y_i$ sont tirés identiquement, indépendamment et aléatoirement de $x_1,\ldots,x_n$ -- aléatoirement signifie ici que chacun de $x_i$ a une probabilité 1/n d'être tiré. Dans l'échantillon de bootstrap $y_1,\ldots,y_n$ nous avons donc certains $x_i$ qui sont présents plusieurs fois alors que d'autres sont absents. Une fois l'échantillon de bootstrap obtenu, la valeur de la statistique d'intérêt est calculée sur celui-ci. On répète ce ré-échantillonage 500 ou 1000 fois et on obtient ainsi 500 ou 1000 valeurs de la statistique d'intérêt. Ces valeurs peuvent être utilisées pour construire un histogramme qui, s'il est correctement normalisé -- sont intégrale est 1 -- servira d'estimation de la loi d'échantillonage de la statistique. Je sais que la première fois qu'on entend parler de cette méthode on se dit que c'est une « escroquerie » mais on peut prouver que ça marche !

La fonction Python qui vous permet de générer un échantillon de bootstrap s'appelle choice et elle se trouve dans le sous-module numpy.random. Vous ferez donc la chose suivante: Pour chacun des 1000 échantillons de bootstrap que vous générerez, vous calculerez $\tilde{a}$, $\widetilde{\mathrm{scale}}$, $\hat{a}$ et $\widehat{\mathrm{scale}}$ -- en utilisant les fonctions que vous aurez définit précédamment -- et vous ferez un histogramme à partir des 1000 valeurs obtenues pour chacune de ces quatre statistiques. Vous ajouterez à vos histogrammes une ligne verticale marquant la valeur de la statistique obtenue sur le vrai échantillon, c'est-à-dire aux questions 2 et 3. Vous utiliserez les percentiles à 2,5 et 97,5 % de ces échantillons de bootstrap pour obtenir des intervalles de confiance à 95 % pour chacune de vos statistiques. Vous n'oublierez pas de spécifier la graine de votre générateur de nombres pseudo-aléatoires avant de lancer vos tirage de bootstrap ; je pourrez ainsi reproduire exactement ce que vous avez fait.

Réponse

Ici, ce qui prend du temps, c'est de lire l'énoncer et de le comprendre ; le code est tout bête :

In [31]:
from numpy.random import choice
np.random.seed(20150618)
nb_rep = 1000
n = len(baleines)
boot_mm = zeros((nb_rep,2))
boot_mv = zeros((nb_rep,3))
for r_idx in range(nb_rep):
    boot_rep = choice(baleines,n)
    boot_mm[r_idx,:] = methode_moments(boot_rep)
    boot_mv[r_idx,:] = methode_mv(boot_rep)

Nous vérifions que notre méthode de Newton a toujours converger avec :

In [32]:
max(boot_mv[:,2])
Out[32]:
2.0

Nous pouvons maintenant faire les histogrammes demandés :

In [37]:
a_tilde_hist = hist(boot_mm[:,0],bins=50,color='black',normed=True)
xlabel(r'$\tilde{a}$')
ylabel(u'Densité estimée')
title(r"Loi d'échantillonage estimée pour $\tilde{a}$")
axvline(a_tilde,color='red',lw=2)
Out[37]:
<matplotlib.lines.Line2D at 0x7fd23d730690>
In [35]:
scale_tilde_hist = hist(boot_mm[:,1],bins=50,color='black',normed=True)
xlabel(r'$\widetilde{\mathrm{scale}}$')
ylabel(u'Densité estimée')
title(r"Loi d'échantillonage estimée pour $\widetilde{\mathrm{scale}}$")
axvline(scale_tilde,color='red',lw=2)
Out[35]:
<matplotlib.lines.Line2D at 0x7fd23d54d5d0>

On remarque que l'estimateur de la méthode des moments pour scale_tilde a l'air biaisé.

In [38]:
a_hat_hist = hist(boot_mv[:,0],bins=50,color='black',normed=True)
xlabel(r'$\hat{a}$')
ylabel(u'Densité estimée')
title(r"Loi d'échantillonage estimée pour $\hat{a}$")
axvline(a_hat,color='red',lw=2)
Out[38]:
<matplotlib.lines.Line2D at 0x7fd23d006490>
In [39]:
scale_hat_hist = hist(boot_mv[:,1],bins=50,color='black',normed=True)
xlabel(r'$\widehat{\mathrm{scale}}$')
ylabel(u'Densité estimée')
title(r"Loi d'échantillonage estimée pour $\widehat{\mathrm{scale}}$")
axvline(scale_hat,color='red',lw=2)
Out[39]:
<matplotlib.lines.Line2D at 0x7fd23cdde050>

Les intervalles de confiance à 95 % sont alors, pour a_tilde :

In [40]:
sort(copy(boot_mm[:,0]))[[25,975]]
Out[40]:
array([ 0.49571648,  0.99503969])

Pour scale_tilde :

In [41]:
sort(copy(boot_mm[:,1]))[[25,975]]
Out[41]:
array([ 0.65092811,  1.1133052 ])

Pour a_hat :

In [42]:
sort(copy(boot_mv[:,0]))[[25,975]]
Out[42]:
array([ 1.35530828,  2.0046141 ])

Pour scale_hat:

In [44]:
sort(copy(boot_mv[:,1]))[[25,975]]
Out[44]:
array([ 0.26250012,  0.51119388])

Test d'adéquation (Question 6, plus difficile -- poids 4 ou 5)

Vous testerez visuellement l'adéquation de la loi gamma avec les paramètres obtenus par la méthode du maximum de vraisemblance en considérant l'histogramme non normalisé -- c'est-à-dire dont chaque bin / classe contient le nombre d'observations -- comme une suite d'observations de variables de Poisson dont le paramètre est donné par la valeur de l'intégrale de la densité de la loi gamma entre les limites de chaque bin, multiplié par le nombre total d'observations. Vous appliquerez alors une stabilisation de variance qui vous permettra de juger graphiquement de l'adéquation de votre loi. Plus explicitement, vous transformerez vos comptes par bin -- suivant la transformation stabilisant la variance --, vous transformerez de même vos prédictions de nombres de comptes par bin obtenues avec votre loi gamma et vous soustrairez la seconde de la première. L'écart de cette différence par-rapport à 0, compte tenu du niveau auquel vous aurez stabilisé la variance, vous indiquera si votre loi gamma ajustée est satisfaisante, ou pas.

Réponse

Ce que je vous demande de faire est en fait un hanging rootogram. Vous trouverez une discussion de ce type d'objet dans mon cours d'il y a deux ans.

Comme nous devons travailler avec les comptes par bin / classe nous devons lire un peu la documentation de la fonction hist pour coprendre ce qu'elle renvoie en plus de son « effet de bord » : la génération du graphe. Là nous voyons qu'elle nous renvoie un vecteur de comptes (si nous n'avons pas normalisé notre histogramme) et un vecteur de frontière entre les bins :

In [45]:
baleine_hist2 = hist(baleines, bins=50, color='black')
In [46]:
len(baleine_hist2)
Out[46]:
3
In [47]:
baleine_hist2[0]
Out[47]:
array([ 33.,  42.,  47.,  22.,   9.,   9.,   4.,   6.,   2.,   2.,   2.,
         3.,   4.,   2.,   2.,   1.,   4.,   0.,   2.,   3.,   0.,   0.,
         1.,   0.,   1.,   1.,   0.,   0.,   2.,   1.,   0.,   0.,   1.,
         0.,   0.,   0.,   1.,   0.,   0.,   0.,   0.,   0.,   1.,   0.,
         0.,   1.,   0.,   0.,   0.,   1.])

Ici nous allons travailler avec la transformation : $2 \times \sqrt{\quad}$. Nous multiplions par 2 afin de stabiliser la variance à 1. D'où :

In [48]:
comptes_stab = 2*sqrt(baleine_hist2[0])

Maintenant l'intégrale de la densité obtenue par EMV entre les frontières des bins n'est rien d'autre que la différence de la fonction de répartition entre ces mêmes frontières soit (attention à ne pas oublier de multiplier par le nombre total d'observations pour avoir la prédiction des comptes):

In [49]:
comptes_pred = diff(gamma.cdf(baleine_hist2[1],a=a_hat,scale=scale_hat))*len(baleines)
In [50]:
comptes_pred_stab = 2*sqrt(comptes_pred)

Nous n'avons plus qu'à faire le graphe :

In [60]:
xx = baleine_hist2[1][:-1]+diff(baleine_hist2[1])/2
plot(xx,comptes_stab-comptes_pred_stab,color='red',lw=2)
xlabel('Temps (heure)')
ylabel(u'Différence obs. pred. après stabilisation de la var. à 1')
axhline(1.96,color='black',linestyle='--')
axhline(-1.96,color='black',linestyle='--')
grid()

Sous l'hypothèse que notre loi ajustée a bien généré les données, on s'attend à ce que sur 50 bins, en moyenne 2,5 soit à l'extérieur de la bandes [-1.96,196] sur le graphe ci-dessus. Ici nous en avons plus que le double. Je me demande vraiment comment des gens ont réussi à publier ce genre d'analyse…

In [ ]: