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

Notebook originel: exam24.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

The history saving thread hit an unexpected error (OperationalError('attempt to write a readonly database')).History will not be written to the database.

[2]:

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

[3]:

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

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🔗

\[\text{MAD} = \text{med}(| x - \text{med}(x) |)\]
\[\hat{\sigma} = \frac{1}{\Phi^{-1}(3/4)}\text{MAD}\]

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.

\[\text{IQR} = \{x\}_{75\%} - \{x\}_{25\%}\]
\[\hat{\sigma} = \frac{1}{2\Phi^{-1}(3/4)} \text{IQR}\]

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):

\[\hat{\sigma} = S_n = c\;\text{med}_i \left( \text{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)\]

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

  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.

  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).

  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).

  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.

  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).

  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).

[5]:

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]])

  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).

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)🔗

[7]:

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_01_std_shape(self):
        with self.assertRaises(AssertionError):
            sigma_std_1D(self.arr2D)       # Invalid shape

    def test_02_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_03_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_04_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_05_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_06_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_07_objfn(self):
        self.assertAlmostEqual(objfn(1), 0.0805853432046475, places=self.places)

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

    def test_09_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_10_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_upper_pairs(self):
    #     np.testing.assert_array_equal(upper_pairs_absdiff(np.arange(3)), [1, 2, 1])

    def test_11_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_12_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_13_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.main(argv=[''], verbosity=2, exit=False);

test_01_std_shape (__main__.TestNotebook) ... ERROR
test_02_std (__main__.TestNotebook) ... ERROR
test_03_mad (__main__.TestNotebook) ... ERROR
test_04_sigma_mad (__main__.TestNotebook) ... ERROR
test_05_iqr (__main__.TestNotebook) ... ERROR
test_06_sigma_iqr (__main__.TestNotebook) ... ERROR
test_07_objfn (__main__.TestNotebook) ... ERROR
test_08_Sn_norm (__main__.TestNotebook) ... ERROR
test_09_all_pairs (__main__.TestNotebook) ... ERROR
test_10_Sn (__main__.TestNotebook) ... ERROR
test_11_Qn (__main__.TestNotebook) ... ERROR
test_12_sigma_mad_nd (__main__.TestNotebook) ... ERROR
test_13_sigma_iqr_nd (__main__.TestNotebook) ... ERROR

======================================================================
ERROR: test_01_std_shape (__main__.TestNotebook)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/tmp/ipykernel_3190279/3804737905.py", line 20, in test_01_std_shape
    sigma_std_1D(self.arr2D)       # Invalid shape
NameError: name 'sigma_std_1D' is not defined

======================================================================
ERROR: test_02_std (__main__.TestNotebook)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/tmp/ipykernel_3190279/3804737905.py", line 23, in test_02_std
    self.assertAlmostEqual(sigma_std_1D(self.arr_no), 0.9697001436773827, places=self.places)
NameError: name 'sigma_std_1D' is not defined

======================================================================
ERROR: test_03_mad (__main__.TestNotebook)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/tmp/ipykernel_3190279/3804737905.py", line 27, in test_03_mad
    self.assertAlmostEqual(mad_1D(self.arr_no), 0.7009412573844651, places=self.places)
NameError: name 'mad_1D' is not defined

======================================================================
ERROR: test_04_sigma_mad (__main__.TestNotebook)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/tmp/ipykernel_3190279/3804737905.py", line 31, in test_04_sigma_mad
    self.assertAlmostEqual(sigma_mad_1D(self.arr_no), 1.0392170632403142, places=self.places)
NameError: name 'sigma_mad_1D' is not defined

======================================================================
ERROR: test_05_iqr (__main__.TestNotebook)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/tmp/ipykernel_3190279/3804737905.py", line 35, in test_05_iqr
    self.assertAlmostEqual(iqr_1D(self.arr_no), 1.4012263172056918, places=self.places)
NameError: name 'iqr_1D' is not defined

======================================================================
ERROR: test_06_sigma_iqr (__main__.TestNotebook)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/tmp/ipykernel_3190279/3804737905.py", line 39, in test_06_sigma_iqr
    self.assertAlmostEqual(sigma_iqr_1D(self.arr_no), 1.0387306232587965, places=self.places)
NameError: name 'sigma_iqr_1D' is not defined

======================================================================
ERROR: test_07_objfn (__main__.TestNotebook)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/tmp/ipykernel_3190279/3804737905.py", line 43, in test_07_objfn
    self.assertAlmostEqual(objfn(1), 0.0805853432046475, places=self.places)
NameError: name 'objfn' is not defined

======================================================================
ERROR: test_08_Sn_norm (__main__.TestNotebook)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/tmp/ipykernel_3190279/3804737905.py", line 46, in test_08_Sn_norm
    self.assertAlmostEqual(Sn_norm(), 1.1925985531232088, places=12)
NameError: name 'Sn_norm' is not defined

======================================================================
ERROR: test_09_all_pairs (__main__.TestNotebook)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/tmp/ipykernel_3190279/3804737905.py", line 49, in test_09_all_pairs
    np.testing.assert_array_equal(all_pairs_absdiff(np.arange(3)), [[0, 1, 2], [1, 0, 1], [2, 1, 0]])
NameError: name 'all_pairs_absdiff' is not defined. Did you mean: 'upper_pairs_absdiff'?

======================================================================
ERROR: test_10_Sn (__main__.TestNotebook)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/tmp/ipykernel_3190279/3804737905.py", line 52, in test_10_Sn
    self.assertAlmostEqual(sigma_Sn_1D(self.arr_no), 0.9833810779051267, places=self.places)
NameError: name 'sigma_Sn_1D' is not defined

======================================================================
ERROR: test_11_Qn (__main__.TestNotebook)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/tmp/ipykernel_3190279/3804737905.py", line 59, in test_11_Qn
    self.assertAlmostEqual(sigma_Qn_1D(self.arr_no), 1.059319495196456, places=self.places)
NameError: name 'sigma_Qn_1D' is not defined

======================================================================
ERROR: test_12_sigma_mad_nd (__main__.TestNotebook)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/tmp/ipykernel_3190279/3804737905.py", line 63, in test_12_sigma_mad_nd
    self.assertAlmostEqual(sigma_mad(self.arr2D), 1.0392170632403142, places=self.places)
NameError: name 'sigma_mad' is not defined

======================================================================
ERROR: test_13_sigma_iqr_nd (__main__.TestNotebook)
----------------------------------------------------------------------
Traceback (most recent call last):
  File "/tmp/ipykernel_3190279/3804737905.py", line 72, in test_13_sigma_iqr_nd
    self.assertAlmostEqual(sigma_iqr(self.arr2D), 1.0387306232587965, places=self.places)
NameError: name 'sigma_iqr' is not defined

----------------------------------------------------------------------
Ran 13 tests in 0.065s

FAILED (errors=13)

[ ]:

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

Contenu