Examen Science des données: estimateurs robustes de la dispersion

Notebook originel: exam24_corrige.ipynb

Examen Science des données: estimateurs robustes de la dispersion🔗

Date: 5 novembre 2024, 10h00-12h00

Durée: 2 h

L’objectif général de l’exercice est de mettre en place des estimateurs robustes de la dispersion d’un échantillon, c.à.d des estimateurs moins sensibles à la présence éventuelle de points abérrants (i.e. des points ne respectant pas la distribution intrinsèque de la population) dans l’échantillon.

Consignes🔗

Créez un répertoire Exam24_NOM_Prenom/ dans lequel vous travaillerez, et où seront stockés

  1. le notebook jupyter d’analyse (nom_prenom.ipynb, copié à partir de ce notebook énoncé sans l’extension .orig), que vous complèterez progressivement en ajoutant une ou plusieurs cellules en dessous de chaque question;

  2. le package python en développement, selon les instructions ci-dessous.

Vous devez alors mettre ce répertoire sous contrôle git sur un projet dédié sur le serveur https://gitlab.in2p3.fr, et en transmettre l’adresse à y.copin@ipnl.in2p3.fr.

Toutefois, si, par manque de temps ou de connaissance, vous ne pouvez pas créer un dépôt git, envoyez par mail le notebook et/ou le package (sous forme d’un fichier tar ou zip) à y.copin@ipnl.in2p3.fr (cette solution sera pénalisée dans la notation).

Tests🔗

Ce notebook inclut des tests unitaires vous permettant de tester a minima le code du notebook. Utilisez ces tests pour guider votre développement et vérifier vos résultats. Ne modifiez pas (a fortiori, n’effacez pas) ces tests.

Si vous n’avez pas réussi à stocker votre package sous ``git``, envoyer votre notebook à y.copin@ipnl.in2p3.fr.

[1]:

%matplotlib inline

[2]:

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats

[3]:

np.set_printoptions(precision=3, suppress=True);

CORRECTION🔗

Référence: https://www.statsmodels.org/dev/_modules/statsmodels/robust/scale.html

Détermination robuste de la dispersion d’un échantillon🔗

Considérons un échantillon monovarié \(x = \{x_i\}_{i=1\ldots N}\). Les statistiques suivantes sont des estimateurs de la dispersion de l’échantillon, renormalisés pour fournir l’écart type dans le cas d’une distribution normale \(\mathcal{N}\).

Écart type🔗

\[\hat{\sigma} = \frac{N}{N - 1} \left<(x - \left<x\right>)^2\right>\]

Median Absolute Deviation🔗

\[\begin{split}\DeclareMathOperator{\med}{med} \begin{align} \text{MAD} &= \med(| x - \med(x) |) \\ \hat{\sigma} &= \frac{1}{\Phi^{-1}(3/4)}\text{MAD} %\simeq 1.4826\times\text{MAD} \end{align}\end{split}\]

où \(\Phi\) est la the fonction de distribution cumulative normale standard (scipy.stats.norm.cdf), et son inverse \(\Phi^{-1}\) est la fonction quantile normale (scipy.stats.norm.ppf).

Écart Inter-Quartile🔗

Il s’agit de la plage tronquée à 25%, c’est-à-dire la différence entre le 75e et le 25e percentile d’un échantillon.

\[\begin{split}\begin{align} \text{IQR} &= \{x\}_{75\%} - \{x\}_{25\%} \\ \hat{\sigma} &= \frac{1}{2\Phi^{-1}(3/4)} \text{IQR} %\simeq 0.7413\times\mathrm{IQR} \end{align}\end{split}\]

Les estimateurs MAD et IQR ont tous deux un point de rupture optimal de 0,5 : jusqu’à la moitié des points de l’échantillon peuvent être arbitrairement grands avant que l’estimateur donne une estimation arbitrairement grande. (En comparaison, l’écart-type de l’échantillon a un point de rupture de 0.) Malheureusement, ces estimateurs sont également moins efficaces statistiquement que l’écart type traditionnel.

Par ailleurs, ces estimateurs sont légèrement biaisés (voir Croux & Rousseeuw, 1992), notamment pour les petits échantillons, mais nous n’aborderons pas ce point dans cet examen.

Différences absolues par pairs🔗

D’autres estimateurs de dispersion avec un point de rupture maximal sont présentés dans Croux and Rousseeuw (1992):

\[\begin{split}\begin{align} \hat{\sigma} = S_n &= c \med_i \left( \med_j(\left| x_i - x_j \right|) \right), \\ \hat{\sigma} = Q_n &= \frac{1}{\sqrt{2}\;\Phi^{-1}(5/8)} \text{first quartile of} \left( \left| x_i - x_j \right| : i < j \right). \end{align}\end{split}\]

Nous les implémenterons dans une version simplifiée.

Échantillons de test (avec et sans outliers)🔗

[4]:

m, n = 10, 8
rng = np.random.default_rng(0)

# Échantillon sans outliers
arr_no = rng.normal(size=m * n)

# Échantillon avec outliers
arr_wo = arr_no.copy()
for i in range(0, m*n, n+1):
    arr_wo[i] += (-1)**i * 6

Questions🔗

  1. Tracer un histogramme normalisé de ces deux échantillons (ax.hist(density=True)), ainsi que la distribution normale \(\mathcal{N}(\mu=0, \sigma=1)\) (scipy.stats.norm.pdf).

    histo.png

[5]:

fig, ax = plt.subplots()
bins = np.linspace(-10, +10, 21)
ax.hist(arr_no, bins=bins, density=True, histtype='step', label="no outliers")
ax.hist(arr_wo, bins=bins, density=True, histtype='step', label="w/ outliers")
xs = np.linspace(-3, +3)
ax.plot(xs, scipy.stats.norm.pdf(xs), label="N(0, 1)")
ax.set(xlabel="Value", ylabel="PDF")
ax.legend()
fig.savefig("histo.png")
../../_images/Annales_24-scale_exam24_corrige_9_0.png

  1. Écrire une fonction sigma_std_1D(arr) qui retourne l’écart type de l’échantillon (np.std(ddof=1)).

    Pour cette fonction et les suivantes, s’assurer (assert) que le tableau d’entrée est unidimensionnel.

[6]:

def sigma_std_1D(arr):

    assert np.ndim(arr) == 1

    return np.std(arr, ddof=1)

[7]:

sigma_std_1D(arr_no), sigma_std_1D(arr_wo),

[7]:

(0.9697001436773827, 2.326360764433049)

  1. Écrire les fonctions mad_1D(arr) et sigma_mad_1D(arr) qui retourne le MAD normalisé de l’échantillon (utiliser np.median et scipy.stats.norm.ppf).

[8]:

def mad_1D(arr):

    assert np.ndim(arr) == 1

    return np.median(np.abs(arr - np.median(arr)))

[9]:

mad_1D(arr_no), mad_1D(arr_wo)

[9]:

(0.7009412573844651, 0.8302102070751392)

[10]:

def sigma_mad_1D(arr, small_sample=False):

    assert np.ndim(arr) == 1

    norm = scipy.stats.norm.ppf(0.75)
    # print(f"{norm=}, {1/norm=}")
    if small_sample:
        # Williams (2011) via Akinshin (2022, 2207.12005)
        bns = [None, None, 1.197, 1.490, 1.360, 1.217, 1.189, 1.138, 1.127, 1.101]
        n, = np.shape(arr)
        if n <= 9:
            bn = bns[n]
        else:
            bn = n / (n - 0.801)
        norm *= bn

    return mad_1D(arr) / norm

[11]:

sigma_mad_1D(arr_no), sigma_mad_1D(arr_wo)

[11]:

(1.0392170632403142, 1.2308714948355965)

  1. Écrire les fonctions iqr_1D(arr) et sigma_iqr_1D(arr) qui retourne l’écart inter-quartile normalisé de l’échantillon (utiliser np.quantile).

[12]:

def iqr_1D(arr):

    assert np.ndim(arr) == 1

    q1, q2 = np.quantile(arr, (.25, .75))

    return q2 - q1

[13]:

iqr_1D(arr_no), iqr_1D(arr_wo)

[13]:

(1.4012263172056918, 1.5901879638726397)

[14]:

def sigma_iqr_1D(arr):

    assert np.ndim(arr) == 1

    norm = 2 * scipy.stats.norm.ppf(0.75)
    # print(f"{norm=}, {1/norm=}")

    return iqr_1D(arr) / norm

[15]:

sigma_iqr_1D(arr_no), sigma_iqr_1D(arr_wo)

[15]:

(1.0387306232587965, 1.1788081015392409)

  1. Pour l’estimateur \(Sn\), il est d’abord nécessaire de déterminer la normalisation \(c\), qui résout cette équation:

    \[\Phi\left(\Phi^{-1}(3/4) + 1/c\right) - \Phi\left(\Phi^{-1}(3/4) - 1/c\right) - 1/2 = 0\]

    1. Définir la fonction objfn(c) dont on cherche le zéro.

    2. Tracer cette fonction sur l’intervale \([1/2, 2]\).

    3. Écrire une fonction Sn_norm() qui retourne la solution de l’équation précédente.

[16]:

def objfn(c):

    return scipy.stats.norm.cdf(scipy.stats.norm.ppf(0.75) + 1/c) - scipy.stats.norm.cdf(scipy.stats.norm.ppf(0.75) - 1/c) - 0.5

[17]:

objfn(1)

[17]:

0.0805853432046475

[18]:

xs = np.geomspace(0.5, 2)
fig, ax = plt.subplots()
ax.plot(xs, objfn(xs))
plt.show();
../../_images/Annales_24-scale_exam24_corrige_26_0.png 
[19]:

def Sn_norm():
    from scipy.optimize import brentq

    zero = brentq(objfn, 0.5, 2)

    return zero

[20]:

Sn_norm()

[20]:

1.1925985531232088

  1. Implémentation (naïve et approximative) de l’estimateur \(Sn\)

    1. Écrire une fonction all_pairs_absdiff(arr) qui retourne le tableau \(N\times N\) des différences absolues \(\left[\,|x_i - x_j|\,\right]_{i,j = 1\ldots N}\) (utiliser une double liste en compréhension).

    2. Écrire la fonction sigma_Sn_1D(arr) qui retourne l’estimateur \(Sn\) normalisé de l’échantillon (utiliser \(c = 1.19259855\) si vous n’avez pas sû répondre à la question 5.C).

[21]:

def all_pairs_absdiff(arr):

    return np.abs([ [ xi - xj for xi in arr ] for xj in arr ])

[22]:

all_pairs_absdiff(np.arange(3))

[22]:

array([[0, 1, 2],
       [1, 0, 1],
       [2, 1, 0]])

[23]:

def sigma_Sn_1D(arr, c=1.19259855):

    assert np.ndim(arr) == 1

    return c * np.median(np.median(all_pairs_absdiff(arr), axis=0))

[24]:

sigma_Sn_1D(arr_no), sigma_Sn_1D(arr_wo)

[24]:

(0.9833810779051267, 1.1717986891987089)

  1. Implémentation (naïve et approximative) de l’estimateur \(Qn\).

    En utilisant la fonction upper_pairs_absdiff(arr) ci-dessous, écrire la fonction sigma_Qn_1D(arr) qui retourne l’estimateur \(Qn\) normalisé de l’échantillon (utiliser np.quantile).

[25]:

def upper_pairs_absdiff(arr):

    n, = np.shape(arr)
    idx = np.triu_indices(n, k=1)

    return np.abs(arr[idx[0]] - arr[idx[1]])

[26]:

upper_pairs_absdiff(np.arange(3))

[26]:

array([1, 2, 1])

[27]:

def sigma_Qn_1D(arr):

    assert np.ndim(arr) == 1

    qn = np.quantile(upper_pairs_absdiff(arr), 0.25)
    norm = 2**0.5 * scipy.stats.norm.ppf(5/8)
    #print(f"{norm=}, {1/norm=}")

    return qn / norm

[28]:

sigma_Qn_1D(arr_no), sigma_Qn_1D(arr_wo)

[28]:

(1.059319495196456, 1.2943757905148814)

  1. Écrire les fonctions sigma_mad(arr, axis=None) et sigma_iqr(arr, axis=None) généralisant les fonctions unidimensionnelles précédentes aux tableaux de rang arbitraire, et en calculant les statistiques le long des axes spécifiés (None pour l’ensemble du tableau).

[29]:

def sigma_mad(arr, axis=None):

    return np.median(np.abs(arr - np.median(arr, axis=axis, keepdims=True)), axis=axis) / scipy.stats.norm.ppf(0.75)

[30]:

sigma_mad(arr_no.reshape(m, n)), sigma_mad(arr_no.reshape(m, n), axis=0), sigma_mad(arr_no.reshape(m, n), axis=1)

[30]:

(1.0392170632403142,
 array([1.073, 0.859, 1.46 , 1.263, 0.395, 0.719, 1.715, 1.262]),
 array([0.573, 0.761, 0.803, 0.727, 0.637, 0.926, 0.978, 1.054, 0.866,
        1.678]))

[31]:

def sigma_iqr(arr, axis=None):

    q1, q2 = np.quantile(arr, (0.25, 0.75), axis=axis)

    return (q2 - q1) / (2 * scipy.stats.norm.ppf(0.75))

[32]:

sigma_iqr(arr_no.reshape(m, n)), sigma_iqr(arr_no.reshape(m, n), axis=0), sigma_iqr(arr_no.reshape(m, n), axis=1)

[32]:

(1.0387306232587965,
 array([0.914, 0.711, 1.174, 1.135, 0.436, 0.647, 1.555, 1.057]),
 array([0.498, 0.54 , 0.699, 0.677, 0.548, 1.002, 0.835, 0.953, 0.644,
        1.362]))

Packaging🔗

  1. Extraire toutes les fonctions précédentes et les inclure dans un fichier scale.py. Il constituera le module de votre package.

  2. Documenter a minima le module et les fonctions par des docstrings (vous pouvez vous aider de l’énoncé).

  3. Créer un package Exam24_Nom_Prenom constitué du seul module scale.py, avec l’arborescence suivante:

    Exam24_Nom_Prenom/
├── scale_nom_prenom
│   ├── scale.py
│   └── __init__.py
├── LICENSE
├── README
├── setup.cfg
└── setup.py

  4. Ajouter un répertoire doc/ pour configurer et générer la documentation Sphinx du package (voir documentation pyyc). Y inclure votre notebook (sous doc/notebooks/) dans une section dédiée.

  5. Ajouter un répertoire test/ pour inclure les tests unitaires de ce notebook (voir documentation pyyc).

  6. Commit/push toutes vos modifications sur votre répo central git, et envoyer l’adresse de ce répo à y.copin@ipnl.in2p3.fr.

    Si vous n’avez pas réussi à stocker votre package sous git, générer une distribution des sources de votre package avec la commande

    $ python setup.py sdist

    et envoyer l’archive correspondante (générée sous dist/) à y.copin@ipnl.in2p3.fr. Si vous n’avez pas non plus réussi à configurer votre package, générer manuellement une archive de votre package:

    $ tar cvzf Exam23_NOM_Prenom.tgz Exam24_NOM_Prenom/

    et envoyer le tarball.

Tests (ne pas modifier)🔗

[35]:

import unittest

class TestNotebook(unittest.TestCase):

    def setUp(self):
        m, n = self.m, self.n = 10, 8
        rng = np.random.default_rng(0)

        self.arr_no = rng.normal(size=m * n)  # Échantillon sans outliers
        self.arr_wo = self.arr_no.copy()      # Échantillon avec outliers
        for i in range(0, m*n, n+1):
            self.arr_wo[i] += (-1)**i * 6

        self.arr2D = self.arr_no.reshape((m, n))

        self.places = 6

    def test_01a_std_shape(self):
        with self.assertRaises(AssertionError):
            sigma_std_1D(self.arr2D)       # Invalid shape

    def test_01b_std(self):
        self.assertAlmostEqual(sigma_std_1D(self.arr_no), 0.9697001436773827, places=self.places)
        self.assertAlmostEqual(sigma_std_1D(self.arr_wo), 2.326360764433049, places=self.places)

    def test_02a_mad_shape(self):
        with self.assertRaises(AssertionError):
            mad_1D(self.arr2D)       # Invalid shape

    def test_02b_mad(self):
        self.assertAlmostEqual(mad_1D(self.arr_no), 0.7009412573844651, places=self.places)
        self.assertAlmostEqual(mad_1D(self.arr_wo), 0.8302102070751392, places=self.places)

    def test_02c_sigma_mad(self):
        self.assertAlmostEqual(sigma_mad_1D(self.arr_no), 1.0392170632403142, places=self.places)
        self.assertAlmostEqual(sigma_mad_1D(self.arr_wo), 1.2308714948355965, places=self.places)

    def test_03a_iqr_shape(self):
        with self.assertRaises(AssertionError):
            iqr_1D(self.arr2D)       # Invalid shape

    def test_03b_iqr(self):
        self.assertAlmostEqual(iqr_1D(self.arr_no), 1.4012263172056918, places=self.places)
        self.assertAlmostEqual(iqr_1D(self.arr_wo), 1.5901879638726397, places=self.places)

    def test_03c_sigma_iqr(self):
        self.assertAlmostEqual(sigma_iqr_1D(self.arr_no), 1.0387306232587965, places=self.places)
        self.assertAlmostEqual(sigma_iqr_1D(self.arr_wo), 1.1788081015392409, places=self.places)

    def test_04a_objfn(self):
        self.assertAlmostEqual(objfn(1), 0.0805853432046475, places=self.places)

    def test_04b_Sn_norm(self):
        self.assertAlmostEqual(Sn_norm(), 1.1925985531232088, places=12)

    def test_05a_all_pairs(self):
        np.testing.assert_array_equal(all_pairs_absdiff(np.arange(3)), [[0, 1, 2], [1, 0, 1], [2, 1, 0]])

    def test_05b_Sn_shape(self):
        with self.assertRaises(AssertionError):
            sigma_Sn_1D(self.arr2D)       # Invalid shape

    def test_05c_Sn(self):
        self.assertAlmostEqual(sigma_Sn_1D(self.arr_no), 0.9833810779051267, places=self.places)
        self.assertAlmostEqual(sigma_Sn_1D(self.arr_wo), 1.1717986891987089, places=self.places)

    def test_06a_upper_pairs(self):
        np.testing.assert_array_equal(upper_pairs_absdiff(np.arange(3)), [1, 2, 1])

    def test_06b_Qn_shape(self):
        with self.assertRaises(AssertionError):
            sigma_Qn_1D(self.arr2D)       # Invalid shape

    def test_06c_Qn(self):
        self.assertAlmostEqual(sigma_Qn_1D(self.arr_no), 1.059319495196456, places=self.places)
        self.assertAlmostEqual(sigma_Qn_1D(self.arr_wo), 1.2943757905148814, places=self.places)

    def test_07a_sigma_mad_nd(self):
        self.assertAlmostEqual(sigma_mad(self.arr2D), 1.0392170632403142, places=self.places)
        np.testing.assert_array_almost_equal(sigma_mad(self.arr2D, axis=0),
                                             [1.073203, 0.858985, 1.459733, 1.263337,
                                              0.395109, 0.719357, 1.715088, 1.262144])
        np.testing.assert_array_almost_equal(sigma_mad(self.arr2D, axis=1),
                                             [0.572676, 0.761405, 0.803225, 0.727297, 0.636682,
                                              0.926232, 0.97819 , 1.053814, 0.866482, 1.677964])

    def test_07b_sigma_iqr_nd(self):
        self.assertAlmostEqual(sigma_iqr(self.arr2D), 1.0387306232587965, places=self.places)
        np.testing.assert_array_almost_equal(sigma_iqr(self.arr2D, axis=0),
                                             [0.914075, 0.710991, 1.173704, 1.135337,
                                              0.436489, 0.647078, 1.554884, 1.056983])
        np.testing.assert_array_almost_equal(sigma_iqr(self.arr2D, axis=1),
                                             [0.497738, 0.540138, 0.698781, 0.677263, 0.547558,
                                              1.002339, 0.835115, 0.952673, 0.643802, 1.361982])

#unittest.TestLoader.sortTestMethodsUsing = lambda *args: +1
unittest.main(argv=[''], verbosity=2, exit=False);

test_01a_std_shape (__main__.TestNotebook) ... ok
test_01b_std (__main__.TestNotebook) ... ok
test_02a_mad_shape (__main__.TestNotebook) ... ok
test_02b_mad (__main__.TestNotebook) ... ok
test_02c_sigma_mad (__main__.TestNotebook) ... ok
test_03a_iqr_shape (__main__.TestNotebook) ... ok
test_03b_iqr (__main__.TestNotebook) ... ok
test_03c_sigma_iqr (__main__.TestNotebook) ... ok
test_04a_objfn (__main__.TestNotebook) ... ok
test_04b_Sn_norm (__main__.TestNotebook) ... ok
test_05a_all_pairs (__main__.TestNotebook) ... ok
test_05b_Sn_shape (__main__.TestNotebook) ... ok
test_05c_Sn (__main__.TestNotebook) ... ok
test_06a_upper_pairs (__main__.TestNotebook) ... ok
test_06b_Qn_shape (__main__.TestNotebook) ... ok
test_06c_Qn (__main__.TestNotebook) ... ok
test_07a_sigma_mad_nd (__main__.TestNotebook) ... ok
test_07b_sigma_iqr_nd (__main__.TestNotebook) ... ok

----------------------------------------------------------------------
Ran 18 tests in 0.022s

OK

[ ]:

Cette page a été générée à partir de exam24_corrige.ipynb.

Contenu