Examen Science des données

Notebook originel: exam22_corrige.ipynb

Examen Science des données🔗

Date: 04 janvier 2023

Durée: 2 h

[1]:

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

Consignes🔗

Créez un répertoire Exam22_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.

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.

[2]:

%matplotlib notebook

Test de Grubbs🔗

En statistique inférentielle, le test de Grubbs (également connu sous le nom de test du résidu normalisé maximum) est un test utilisé pour détecter une seule valeur aberrante dans un ensemble de données univariées supposées provenir d’une population approximativement normale. (NIST/SEMATECH e-Handbook of Statistical Methods).

Attention:, en raison de potentiels effets d’écrantage (plusieurs valeurs aberrantes peuvent affaiblir la performance du test), il n’est généralement pas approprié d’appliquer séquentiellement un test de valeur aberrante unique pour détecter plusieurs valeurs aberrantes.

Le test de Grubbs est défini pour l’hypothèse suivante:

  • \(H_{0}\) : il n’y a pas de valeurs aberrantes dans l’ensemble de données;

  • \(H_{1}\) : il y a exactement une valeur aberrante dans l’ensemble de données.

La statistique du test de Grubbs est définie par :

\[G = \max_{i}\left|\frac{y_{i} − \bar{y}}{s}\right|\]

où \(\bar{y}\) et \(s\) désignent respectivement la moyenne et l’écart type de l’échantillon. Ainsi, \(G\) est le plus grand écart absolu par rapport à la moyenne de l’échantillon en unités de l’écart type de l’échantillon.

Pour le test bilatéral (sans a priori sur la position relative du point abérrant), l’hypothèse \(H_{0}\) est rejetée au niveau de signification \(\alpha\) si

\[G > G_{\max}(N, \alpha) \triangleq \frac{N-1}{\sqrt{N}} \sqrt{\dfrac{\left(t^{*}_{\alpha /(2N),N-2}\right)^{2}}{N - 2 + \left(t^{*}_{\alpha/(2N),N-2}\right)^{2}}}\]

où \(t^{*}_{\alpha/(2N),N-2}\) désigne la valeur critique supérieure de la distribution de Student \(t\) avec \(N - 2\) degrés de liberté et un niveau de signification de \(\alpha/(2N)\).

Dans les cas unilatéraux (tests portant sur \((y_{\max} - \bar{y})/s\) pour la valeur maximale ou \((\bar{y} - y_{\min})/s\) pour la valeur minimale), il faut remplacer \(\alpha/(2N)\) par \(\alpha/N\).

  1. Écrire une fonction zscore(y) agissant sur un tableau \(y\) et retournant le score \(z\) pour chaque élément:

    \[z_{i} = \frac{y_{i} - \bar{y}}{s}\]

    où \(\bar{y}\) et \(s\) s désignent respectivement la moyenne (numpy.mean) et l’écart type (numpy.std) de l’échantillon.

[3]:

def zscore(y, ddof=0):
    """
    Compute the z-score of observations *y* (raveled).

    >>> y = np.array([201.92, 202.18, 200.82, 201.95, 200.19, 245.57, 199.31, 199.53])
    >>> zscore(y)
    array([-0.3043925 , -0.28685895, -0.37857289, -0.3023694 ,
           -0.42105803,  2.63922038, -0.48040234, -0.46556626])
    """

    m = np.mean(y)              # Mean
    s = np.std(y, ddof=ddof)    # StdDev

    return (y.ravel() - m) / s  # z-score

[4]:

y = np.array([201.92, 202.18, 200.82, 201.95, 200.19, 245.57, 199.31, 199.53])
zscore(y)

[4]:

array([-0.3043925 , -0.28685895, -0.37857289, -0.3023694 , -0.42105803,
        2.63922038, -0.48040234, -0.46556626])

  1. Écrire une fonction grubbs_gmax(n, alpha=0.05, onesided=False) retournant la valeur critique \(G_{\max}(N, \alpha)\) de la statistique de Grubbs pour \(N\) points de données et un niveau de signification \(\alpha\), dans le cas bilatéral (onesided=False) ou unilatéral (True).

    Utiliser scipy.stats.t.isf(epsilon, m) pour la valeur critique \(t^*_{\epsilon, M}\) de la distribution de Student: \(\int_{t^*}^{\infty} t_{m}(x) \text{d}x = \epsilon\).

[5]:

def grubbs_gmax(n, alpha=0.05, onesided=False):
    """
    Critical value for Grubbs' test for outlier.

    >>> grubbs_gmax(n=8, alpha=0.05, onesided=False)
    2.126645087195626
    >>> grubbs_gmax(n=8, alpha=0.05, onesided=True)
    2.031652001549952
    """

    sl = alpha / n    # one-sided
    if not onesided:  # two-sided
        sl /= 2

    # Critical value (squared) of the t distribution w/ n-2 DoF and
    # significance level sl
    cv2 = scipy.stats.t.isf(sl, n - 2) ** 2
    gmax = (n - 1) * np.sqrt(cv2 / (n - 2 + cv2) / n)

    return gmax

[6]:

fig, ax = plt.subplots(tight_layout=True)
alphas = np.geomspace(1e-3, 1e-1, 21)
ns = (5, 10, 20, 50, 100)
for n in ns:
    l, = ax.plot(alphas, [ grubbs_gmax(n, alpha=alpha, onesided=False) for alpha in alphas ],
                 label=f"N={n}")
    ax.plot(alphas, [ grubbs_gmax(n, alpha=alpha, onesided=True) for alpha in alphas ],
            c=l.get_color(), ls='--')
ax.legend()
ax.grid(which='both')
ax.set(xlabel=r"$\alpha$", xscale='log',
       ylabel=r"$G_{\max}$",
       title="Grubbs test critical value (—: two-sided, --: one-sided)");

  1. Écrire une fonction grubbs_test_twosided(y, alpha=0.05) qui:

    1. calcule le score \(z_i\) pour chacun des \(N\) points de \(y\),

    2. calcule la valeur critique \(G_{\max}(N, \alpha)\) (cas bilatéral),

    3. calcule la statistique \(G\) (utiliser numpy.amax),

    4. retourne le tableau de booléens (z == g) & (g > gmax) pour chacun des points.

    Un seul de tous les points peut éventuellement passer ce test et être qualifié d’abérrant.

[7]:

def grubbs_test_twosided(y, alpha=0.05, verbose=False):

    n = len(y.ravel())
    gmax = grubbs_gmax(n, alpha, onesided=False)

    absz = np.abs(zscore(y.ravel()))
    # Grubbs' test statistic
    g = np.amax(absz)
    is_outlier = (absz == g) & (absz > gmax)

    return is_outlier

  1. Écrire une fonction grubbs_test(y, side=0, alpha=0.05) généralisant la fonction précédente aux cas unilatéraux (side=+1 pour le test du maximum, -1 pour le minimum, side=0 pour le test bilatéral). Attention aux cas ``side=±1``.

[8]:

def grubbs_test(y, side=0, alpha=0.05, verbose=False):
    """
    Grubbs' test for outliers (Grubbs 1969 and Stefansky 1972, also known as
    *maximum normed residual test*) is used to detect a *single* outlier in a
    univariate data set that follows an approximately normal distribution.

    Source: `NIST/SEMATECH e-Handbook of Statistical Methods
    <http://www.itl.nist.gov/div898/handbook/eda/section3/eda35h.htm>`_

    .. Warning:: Because of potential screening effects (multiple outliers
       weaken the single-outlier test performance), it is not appropriate to
       apply a single-outlier test sequentially in order to detect multiple
       outliers.

    >>> y = np.array([201.92, 202.18, 200.82, 201.95, 200.19, 245.57, 199.31, 199.53])
    >>> grubbs_test(y, side=0, alpha=0.05)
    array([False, False, False, False, False,  True, False, False])
    >>> grubbs_test(y, side=+1, alpha=0.05)
    array([False, False, False, False, False,  True, False, False])
    >>> grubbs_test(y, side=-1, alpha=0.05)
    array([False, False, False, False, False, False, False, False])
    """

    if side not in (-1, 0, +1):
        raise NotImplementedError(
            "Grubbs' test for outliers is implemented for side=0 and +/-1.")

    n = len(y.ravel())
    gmax = grubbs_gmax(n, alpha, onesided=bool(side))

    z = zscore(y.ravel())

    # Grubbs' test statistic
    if side == 0:               # Two-sided test
        absz = np.abs(z)
        g = np.amax(absz)
        adj = 'extremum'
        label = 'a two-tailed'
        is_outlier = (absz == g) & (absz > gmax)
    elif side == -1:            # One-sided test on minimal value
        g = np.amin(z)
        adj = 'minimum'
        label = 'a lower one-tailed'
        is_outlier = (z == g) & (z < -gmax)
    elif side == +1:            # One-sided test on maximal value
        g = np.amax(z)
        adj = 'maximum'
        label = 'an upper one-tailed'
        is_outlier = (z == g) & (z > gmax)

    if verbose:
        print(f"""\
H0:  there are no outliers in the data
H1:  the {adj} value is an outlier

Test statistic:      G = {g}
Significance level:  α = {alpha}
Critical value for {label} test: G_max = {gmax}
Critical region:     Reject H0 if G > G_max""")

    return is_outlier     # Outlier flags (w/ at most *one* identified outlier)


def grubbs_test_twosided(y, alpha=0.05, verbose=False):
    """
    Two-sided version of `grubbs_test`.
    """

    return grubbs_test(y, side=0, alpha=alpha, verbose=verbose)

[9]:

y = np.array([201.92, 202.18, 200.82, 201.95, 200.19, 245.57, 199.31, 199.53])
grubbs_test(-y, side=0, alpha=0.05, verbose=True)

H0:  there are no outliers in the data
H1:  the extremum value is an outlier

Test statistic:      G = 2.6392203839283317
Significance level:  α = 0.05
Critical value for a two-tailed test: G_max = 2.1266450871956275
Critical region:     Reject H0 if G > G_max

[9]:

array([False, False, False, False, False,  True, False, False])

[10]:

grubbs_test(y, side=+1, alpha=0.05, verbose=True)

H0:  there are no outliers in the data
H1:  the maximum value is an outlier

Test statistic:      G = 2.6392203839283317
Significance level:  α = 0.05
Critical value for an upper one-tailed test: G_max = 2.031652001549951
Critical region:     Reject H0 if G > G_max

[10]:

array([False, False, False, False, False,  True, False, False])

[11]:

grubbs_test(y, side=-1, alpha=0.05, verbose=True)

H0:  there are no outliers in the data
H1:  the minimum value is an outlier

Test statistic:      G = -0.48040234335199367
Significance level:  α = 0.05
Critical value for a lower one-tailed test: G_max = 2.031652001549951
Critical region:     Reject H0 if G > G_max

[11]:

array([False, False, False, False, False, False, False, False])

  1. Écrire une fonction grubbs_clip(y, side=0, alpha=0.05) retournant le jeu de données \(y\) expurgé de son éventuel point abérrant d’après le test de Grubbs.

[12]:

def grubbs_clip(y, side=0, alpha=0.05):
    """
    Remove the outlier (if any) as detected by Grubbs' test.

    >>> y = np.array([201.92, 202.18, 200.82, 201.95, 200.19, 245.57, 199.31, 199.53])
    >>> grubbs_clip(y, side=0, alpha=0.05)
    array([201.92, 202.18, 200.82, 201.95, 200.19, 199.31, 199.53])
    """

    outlier = grubbs_test(y, side=side, alpha=alpha)

    return y[~outlier] # or np.delete(y, outlier)

[13]:

y = np.array([201.92, 202.18, 200.82, 201.95, 200.19, 245.57, 199.31, 199.53])
grubbs_clip(y, side=0, alpha=0.05)

[13]:

array([201.92, 202.18, 200.82, 201.95, 200.19, 199.31, 199.53])

Droite de Henry🔗

Le test de Grubbs est basé sur l’hypothèse de normalité: il faut s’assurer que les données (hors points abérrants) sont raisonnablement approximées par une distribution normale pour l’appliquer.

La droite de Henry est une méthode graphique de type « diagramme Quantile-Quantile » pour estimer la normalité d’une série d’observations. Elle compare la distribution observée des quantiles \(o_{i}\) (« ordered response » en ordonnées, càd juste la version ordonnée de \(y\)) à la distribution attendue des quantiles \(z_{i}\) (« order statistics median » en absisses) dans le cas d’une distribution normale, p.ex. (Wikipedia):

\[\begin{split}\begin{aligned} z_{i} &= \Phi^{-1}(q_{i}) \\ q_{i} &= \frac{i - a}{n + 1 - 2a} \quad i=1\ldots N\\ a &= \begin{cases} 3/8 \quad &N \leq 10 \\ 1/2 \quad &N > 10 \end{cases} \end{aligned}\end{split}\]

où \(\Phi^{-1}\) est la fonction quantile normale (scipy.stats.norm.ppf).

Pour information, cette définition des « rankits » \(q\) est légèrement différente de l’implémentation de scipy.stats.probplot.

  1. Écrire une fonction normal_quantiles(y) retournant \((z, o)\) définis ci-dessus pour un jeu de données \(y\).

[14]:

def normal_quantiles(y, rankit='R'):
    """
    Compute the *ordered response* and the *order statistics median* for dataset *y*.

    >>> y = np.array([201.92, 202.18, 200.82, 201.95, 200.19, 245.57, 199.31, 199.53])
    >>> normal_quantiles(y)
    (array([-1.43420016, -0.85249503, -0.47278912, -0.15250597,
             0.15250597,  0.47278912,  0.85249503,  1.43420016]),
     array([199.31, 199.53, 200.19, 200.82, 201.92, 201.95, 202.18, 245.57]))
    """

    n = len(y.ravel())
    i = np.arange(n) + 1
    if rankit == 'R':
        # Définition R via Wikipedia (https://en.wikipedia.org/wiki/Normal_probability_plot#Definition)
        a = 3/8 if n <= 10 else 0.5
        q = (i - a) / (n + 1 - 2*a)
    elif rankit == 'scipy':
        # Définition scipy.stats.probplot
        q = (i - 0.3175) / (n + 0.365)
        q[-1] = 0.5 ** (1/n)
        q[0] = 1 - q[-1]
    else:
        raise NotImplementedError(f"Unknown rankit definition {rankit!r}.")

    osm = scipy.stats.norm.ppf(q)
    osr = np.sort(y)

    return osm, osr

[15]:

y = np.array([201.92, 202.18, 200.82, 201.95, 200.19, 245.57, 199.31, 199.53])
normal_quantiles(y)

[15]:

(array([-1.43420016, -0.85249503, -0.47278912, -0.15250597,  0.15250597,
         0.47278912,  0.85249503,  1.43420016]),
 array([199.31, 199.53, 200.19, 200.82, 201.92, 201.95, 202.18, 245.57]))

  1. Écrire une fonction graphique plot_quantiles(y) traçant \(o\) en fonction de \(z\), et y ajoutant la droite de régression linéaire et le coefficient de corrélation linéaire (utiliser scipy.stats.linregress), p.ex.:

    >>> rng = np.random.default_rng(0)
>>> y = rng.normal(size=100)
>>> plot_quantiles(y)

    henry.png

[16]:

def plot_quantiles(y, ax=None, alpha=0.05):

    osm, osr = normal_quantiles(y)
    result = scipy.stats.linregress(osm, osr)

    c = 'C00' if result.pvalue < alpha else 'C03'  # bleu si normal, rouge sinon

    if ax is None:
        fig, ax = plt.subplots()

    ax.scatter(osm, osr, marker='+', c='k', label="Data")
    ax.plot(osm[[0, -1]], result.intercept + result.slope * osm[[0, -1]],
            label=f"r={result.rvalue:.3f}", c=c)
    ax.set(xlabel="Expected quantile",
           ylabel="Observed quantile")
    ax.legend()

    return ax

[17]:

rng = np.random.default_rng(0)
y = rng.normal(size=100)
osm, osr = normal_quantiles(y)
result = scipy.stats.linregress(osm, osr)
print(result)

LinregressResult(slope=0.9605701859862437, intercept=0.08109669349071556, rvalue=0.9920231919350392, pvalue=5.803370920219443e-90, stderr=0.01232979926297222, intercept_stderr=0.01225131462706626)

[18]:

plot_quantiles(y)

[18]:

<AxesSubplot: xlabel='Expected quantile', ylabel='Observed quantile'>

  1. Illustrer la pertinence du test de Grubbs sur le jeu de données suivant en traçant les droites de Henry (plot_quantiles) avant puis après élimination du point abérrant:

    >>> y = np.array([201.92, 202.18, 200.82, 201.95, 200.19, 245.57, 199.31, 199.53])

[19]:

y = np.array([201.92, 202.18, 200.82, 201.95, 200.19, 245.57, 199.31, 199.53])
y_clean = grubbs_clip(y)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)
plot_quantiles(y, ax=ax1)
ax1.set_title("Before clipping")
plot_quantiles(y_clean, ax=ax2)
ax2.set(title="After clipping", ylabel='');

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

Packaging🔗

  1. Extraire les fonctions précédentes (zscore, grubbs_gmax, grubbs_test_twosided, grubbs_test, grubbs_clip, normal_quantiles et plot_quantiles) et les inclure dans un fichier grubbs.py. Il constituera le module de votre package.

  2. Documenter le module et les fonctions par des docstrings (objectif, paramètres, valeur de retour, etc.).

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

    Exam22_Nom_Prenom/
├── grubbs_nom_prenom
│   ├── grubbs.py
│   └── __init__.py
├── LICENSE
├── README
├── setup.cfg
└── setup.py

    avec les fichiers suivants (à adapter à votre contexte):

  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. Si vous n’avez pas réussi à stocker votre package sous git, générer un tarball des sources de votre package (duquel vous aurez effacé les répertoires _build/, dist/ s’ils existent),

    $ tar cvzf Exam22_NOM_Prenom.tgz Exam22_NOM_Prenom/

    et envoyer le tarball résultant à nouveau à y.copin@ipnl.in2p3.fr.

Tests🔗

Tests unitaires🔗

[20]:

import unittest

class TestNotebook(unittest.TestCase):

    def setUp(self):
        self.y = np.array([201.92, 202.18, 200.82, 201.95, 200.19, 245.57, 199.31, 199.53])

    def test_01_zscore(self):
        np.testing.assert_allclose(
            zscore(self.y),
            [-0.3043925 , -0.28685895, -0.37857289, -0.3023694 ,
             -0.42105803,  2.63922038, -0.48040234, -0.46556626])

    def test_02_grubbs_gmax(self):
        self.assertAlmostEqual(grubbs_gmax(n=8, alpha=0.05, onesided=False),
                               2.126645087195626)
        self.assertAlmostEqual(grubbs_gmax(n=8, alpha=0.05, onesided=True),
                               2.031652001549952)
        self.assertAlmostEqual(grubbs_gmax(n=8, alpha=0.01, onesided=False),
                               2.2743651271139984)

    def test_03_grubbs_test_twosided(self):
        np.testing.assert_equal(grubbs_test_twosided(self.y, alpha=0.05),
                                [False, False, False, False, False,  True, False, False])
        np.testing.assert_equal(grubbs_test_twosided(-self.y, alpha=0.05),
                                [False, False, False, False, False,  True, False, False])

    def test_04a_grubbs_test_both(self):
        np.testing.assert_equal(grubbs_test(self.y, side=0, alpha=0.05),
                                [False, False, False, False, False,  True, False, False])
        np.testing.assert_equal(grubbs_test(-self.y, side=0, alpha=0.05),
                                [False, False, False, False, False,  True, False, False])

    def test_04b_grubbs_test_high(self):
        np.testing.assert_equal(grubbs_test(self.y, side=+1, alpha=0.05),
                                [False, False, False, False, False,  True, False, False])
        np.testing.assert_equal(grubbs_test(-self.y, side=+1, alpha=0.05),
                                [False, False, False, False, False, False, False, False])

    def test_04c_grubbs_test_low(self):
        np.testing.assert_equal(grubbs_test(self.y, side=-1, alpha=0.05),
                                [False, False, False, False, False, False, False, False])
        np.testing.assert_equal(grubbs_test(-self.y, side=-1, alpha=0.05),
                                [False, False, False, False, False,  True, False, False])

    def  test_05_grubbs_clip(self):
        np.testing.assert_allclose(
            grubbs_clip(self.y, side=0, alpha=0.05),
            [201.92, 202.18, 200.82, 201.95, 200.19, 199.31, 199.53])

    def test_06a_normal_quantiles_osm_lown(self):
        osm, osr = normal_quantiles(self.y)                # 10 elements < 10
        np.testing.assert_allclose(
            osm,
            [-1.43420016, -0.85249503, -0.47278912, -0.15250597,
              0.15250597,  0.47278912,  0.85249503,  1.43420016])

    def test_06b_normal_quantiles_osm_highn(self):
        y = np.concatenate((self.y, np.diff(self.y[:5])))  # 12 elements > 10
        osm, osr = normal_quantiles(y)
        np.testing.assert_allclose(
            osm,
            [-1.7316644, -1.15034938, -0.8122178, -0.54852228, -0.31863936, -0.10463346,
              0.10463346, 0.31863936,  0.54852228, 0.8122178,   1.15034938,  1.7316644 ])

    def test_06c_normal_quantiles_osr(self):
        osm, osr = normal_quantiles(self.y)
        np.testing.assert_allclose(
            osr,
            [199.31, 199.53, 200.19, 200.82, 201.92, 201.95, 202.18, 245.57])

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

test_01_zscore (__main__.TestNotebook) ... ok
test_02_grubbs_gmax (__main__.TestNotebook) ... ok
test_03_grubbs_test_twosided (__main__.TestNotebook) ... ok
test_04a_grubbs_test_both (__main__.TestNotebook) ... ok
test_04b_grubbs_test_high (__main__.TestNotebook) ... ok
test_04c_grubbs_test_low (__main__.TestNotebook) ... ok
test_05_grubbs_clip (__main__.TestNotebook) ... ok
test_06a_normal_quantiles_osm_lown (__main__.TestNotebook) ... ok
test_06b_normal_quantiles_osm_highn (__main__.TestNotebook) ... ok
test_06c_normal_quantiles_osr (__main__.TestNotebook) ... ok

----------------------------------------------------------------------
Ran 10 tests in 0.011s

OK

[20]:

<unittest.main.TestProgram at 0x7f5431bfb040>

[ ]:

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

Contenu