Notebook originel: exam22_corrige.ipynb

Examen Science des données

Date: 08 décembre 2022

Durée: 2 h

import numpy as np
import matplotlib.pyplot as plt

%matplotlib notebook

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.

Régression linéaire

Soit un jeu de données de \(N\) points \(\{(x_i, y_i)\}\). La régression linéaire

\[y \equiv a + b x\]

est définie par

\[\begin{split}\begin{aligned} \hat{b} &= \frac{\mathrm{Cov}(x, y)}{\mathrm{Var}(x)} \\ \hat{a} &= \langle y\rangle - \hat{b}\langle x\rangle \end{aligned}\end{split}\]

où

\[\begin{split}\begin{aligned} \mathrm{Var}(x) &= \langle(x - \langle x\rangle)^2\rangle \\ \mathrm{Cov}(x, y) &= \langle(x - \langle x\rangle)(y - \langle y\rangle)\rangle \end{aligned}\end{split}\]

Le coefficient de corrélation linéaire est défini par:

\[r = \frac{\mathrm{Cov}(x, y)}{\sqrt{\mathrm{Var}(x)\mathrm{Var}(y)}}\]

1. Sans détailler les calculs, quelles équations ont conduit à ces estimateurs?

2. Dans le notebook, écrire une fonction linregress_1D(x, y) qui, en utilisant uniquement la fonction np.mean, retourne \(\hat{a}\), \(\hat{b}\) et \(r\) pour des vecteurs (1D) d’entrée \(x\) et \(y\).

def linregress_1D(x, y):
    """
    Régression linéaire `y = a + b x`.

    :param np.ndarray x: vecteur des variables explicatives
    :param np.ndarray y: vecteur des variables expliquées
    :return: `a` (ordonnée à l'origine), `b` (pente), `r` (coefficient de corrélation)

    >>> x = np.arange(10)
    >>> linregress_1D(x, 3 + 2 * x)
    (3.0, 2.0, 1.0)
    >>> linregress_1D(x, 3 + 2 * x + np.round(np.sin(x * 10), 3))
    (2.8803090909090923, 2.0317757575757573, 0.9925669351010524)
    """

    xbar = np.mean(x)
    ybar = np.mean(y)
    varx = np.mean((x - xbar)**2)
    vary = np.mean((y - ybar)**2)
    covxy = np.mean((x - xbar)*(y - ybar))

    b = covxy / varx
    a = ybar - b * xbar
    r = covxy / (varx * vary)**0.5

    return a, b, r

3. Tester explicitement sur le cas

>>> x = np.arange(10)
>>> y = 3. + 2. * x + np.round(np.sin(x * 10), 3)

et vérifier le résultat avec scipy.stats.linregress.

n, a, b = 10, 3., 2.
x = np.arange(n)
y = a + b * x
y += np.round(np.sin(x * 10), 3)

res = aa, ba, ra = linregress_1D(x, y)
print(res)

(2.8803090909090923, 2.0317757575757573, 0.9925669351010524)

from scipy.stats import linregress

res = linregress(x, y)
print(res)

LinregressResult(slope=2.031775757575758, intercept=2.8803090909090905, rvalue=0.9925669351010524, pvalue=1.323642496012334e-08, stderr=0.08807677394385144, intercept_stderr=0.47020130379248715)

4. Faire une figure illustrant les points de données \(\{(x, y)\}\) (scatter) et la régression linéaire (plot). Ajouter un label sur chacun des axes, un titre à la figure et une légende.

fig, ax = plt.subplots()
ax.scatter(x, y, c='k')
ax.plot(x, a + b * x, label=f"Input: a={a}, b={b}")
ax.plot(x, aa + ba * x, label=f"Analytic: a={aa:.2f}, b={ba:.2f} (r={ra:.3f})")
ax.set(xlabel='x', ylabel='y')
ax.legend()

<matplotlib.legend.Legend at 0x7f7a291bf130>

Detrending

Le detrending permet de retirer une composant linéaire d’un jeu de données \(\{y_i\}\):

\[\tilde{y}_i = y_i - (\hat{a} + \hat{b} i)\]

1. Écrire une fonction detrend_1D(y) retournant \(\tilde{y}\) (càd \(y\) duquel a été soustrait sa composante linéaire) et \(r\).

def detrend_1D(y):
    """
    Soustraction de la composante linéaire du vecteur (1D) `y`:

    .. math::

       \tilde{y}_i = y_i - (\hat{a} + \hat{b} i)

    :param np.ndarray y: vecteur d'entrée
    :return: $\tilde{y}

    >>> x = np.arange(5)
    >>> detrend_1D(1 + 2 * x)
    array([0., 0., 0., 0., 0.])
    >>> detrend_1D(x**0.5)
    array([-0.28284271,  0.24395221,  0.18496069,  0.02959285, -0.17566304])
    """

    i = np.arange(len(y))
    a, b, r = linregress_1D(i, y)

    return y - (a + b * i)

2. Tester et faire une figure sur le vecteur y précédent.

ytilde = detrend_1D(y)

fig, ax = plt.subplots()
ax.scatter(range(len(y)), ytilde, c='k')
ax.plot(0 * y)
ax.set(xlabel='i', ylabel=r'$\tilde{y}$');

Généralisation multi-dimensionnelle

1. Écrire la function de régression linéaire linregress(x, y, axis) pour retourner le résultat de la régression sur les tableaux x et y de rang arbitraire, le long de la dimension axis. La fonction doit retourner 3 tableaux a, b et r de même rang que les tableaux d’entrée.

   Par example:

   >>> x = np.arange(3*4).reshape((3, 4))
   >>> y = 3 + 2 * x
   >>> linregress(x, y, axis=0)
   (array([[3., 3., 3., 3.]]),
    array([[2., 2., 2., 2.]]),
    array([[1., 1., 1., 1.]]))
   >>> linregress(x, y, axis=1)
   (array([[3.],
           [3.],
           [3.]]),
    array([[2.],
           [2.],
           [2.]]),
    array([[1.],
           [1.],
           [1.]]))
   >>> linregress(x, y, axis=2)
   Traceback (most recent call last):
       ...
   AxisError: axis 2 is out of bounds for array of dimension 2
   

def linregress(x, y, axis=None):
    """
    Régression linéaire `y = a + b x` le long de la dimension `axis`.

    :param np.ndarray x: vecteur des variables explicatives
    :param np.ndarray y: vecteur des variables expliquées
    :param int axis: dimension le long de laquelle réaliser la régression linéaire
    :return: `a` (ordonnée à l'origine), `b` (pente), `r` (coefficient de corrélation)

    >>> x = np.arange(3*4).reshape((3, 4))
    >>> y = 3 + 2 * x
    >>> linregress(x, y, axis=0)  # doctest: +NORMALIZE_WHITESPACE
    (array([[3., 3., 3., 3.]]),
     array([[2., 2., 2., 2.]]),
     array([[1., 1., 1., 1.]]))
    >>> linregress(x, y, axis=1)  # doctest: +NORMALIZE_WHITESPACE
    (array([[3.], [3.], [3.]]),
     array([[2.], [2.], [2.]]),
     array([[1.], [1.], [1.]]))
    >>> linregress(x, y, axis=2)
    Traceback (most recent call last):
        ...
    numpy.AxisError: axis 2 is out of bounds for array of dimension 2
    """

    if axis is None:
        return linregress_1D(np.ravel(x), np.ravel(y))

    xbar = np.mean(x, axis=axis, keepdims=True)
    ybar = np.mean(y, axis=axis, keepdims=True)
    varx = np.mean((x - xbar)**2, axis=axis, keepdims=True)
    vary = np.mean((y - ybar)**2, axis=axis, keepdims=True)
    covxy = np.mean((x - xbar)*(y - ybar), axis=axis, keepdims=True)

    b = covxy / varx
    a = ybar - b * xbar
    r = covxy / (varx * vary)**0.5

    return a, b, r

x = np.arange(3*4).reshape((3, 4))
y = 3 + 2 * x
linregress(x, y, axis=0)

(array([[3., 3., 3., 3.]]),
 array([[2., 2., 2., 2.]]),
 array([[1., 1., 1., 1.]]))

linregress(x, y, axis=1)

(array([[3.],
        [3.],
        [3.]]),
 array([[2.],
        [2.],
        [2.]]),
 array([[1.],
        [1.],
        [1.]]))

2. De même, écrire la fonction detrend(y, axis). Vous pouvez utiliser les lignes suivantes pour générer un indice courant le long de la dimension axis:

   shape = np.ones_like(np.shape(y))
   shape[axis] = -1
   i = np.arange(np.shape(y)[axis]).reshape(shape)
   

def detrend(y, axis=None):
    """
    Soustraction de la composante linéaire `y - (a + b i)` le long de la dimension `axis`.

    :param np.ndarray y: tableau d'entrée
    :param int axis: dimension le long de laquelle soustraire la composante linéaire
    :return: le tableau sans sa composante linéaire

    >>> x = np.arange(3*4).reshape((3, 4))
    >>> y = 3 + 2 * x
    >>> detrend(y, axis=0)  # doctest: +NORMALIZE_WHITESPACE
    array([[0., 0., 0., 0.],
           [0., 0., 0., 0.],
           [0., 0., 0., 0.]])
    >>> detrend(y, axis=1)  # doctest: +NORMALIZE_WHITESPACE
    array([[0., 0., 0., 0.],
           [0., 0., 0., 0.],
           [0., 0., 0., 0.]])
    >>> detrend(y, axis=2)
    Traceback (most recent call last):
        ...
    IndexError: index 2 is out of bounds for axis 0 with size 2
    """

    if axis is None:
        return detrend_1D(np.array(y).ravel()).reshape(np.shape(y))

    shape = np.ones_like(np.shape(y))
    shape[axis] = -1
    i = np.arange(np.shape(y)[axis]).reshape(shape)
    a, b, r = linregress(i, y, axis=axis)

    return y - (a + b * i)

x = np.arange(3*4).reshape((3, 4))
y = np.arange(3).reshape(-1, 1) + 2 * np.arange(1, 4).reshape(-1, 1) * x
y = 3 + 2 * x
print(y, '\n', detrend(y, axis=0), sep='')

[[ 3  5  7  9]
 [11 13 15 17]
 [19 21 23 25]]
[[0. 0. 0. 0.]
 [0. 0. 0. 0.]
 [0. 0. 0. 0.]]

3. Ajouter des options suivantes à la fonction detrend (voir scipy.signal.detrend pour un modèle):

   • type='linear' pour soustraire la composante linéaire \(a + b\,i\) (type='linear'), la composante affine \(b\,i\) (type='affine'), ou uniquement la composante constante \(a\) (type='constant');

   • overwrite=False pour réaliser la soustraction sur un nouveau tableau (overwrite=False) ou sur place (overwrite=True)

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

Packaging

1. Extraire les 4 fonctions précédentes (linregress_1D, detrend_1D, linregress et detrend) et les inclure dans un fichier detrend.py. Il constituera le module de votre package.

2. Documenter (succinctement) le module et les fonctions par des docstrings.

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

   Exam22_Nom_Prenom/
   ├── detrend_nom_prenom
   │   ├── detrend.py
   │   └── __init__.py
   ├── LICENSE
   ├── README
   ├── setup.cfg
   └── setup.py
   

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

   • setup.cfg

   • setup.py

   • __init__.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. 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

import unittest

class TestNotebook(unittest.TestCase):

    def test_linregress_1D(self):
        x = np.arange(10)
        self.assertEqual(linregress_1D(x, 3 + 2 * x), (3.0, 2.0, 1.0))
        self.assertAlmostEqual(linregress_1D(x, 3 + 2 * x + np.round(np.sin(x * 10), 3)),
                               (2.8803090909090923, 2.0317757575757573, 0.9925669351010524))

    def test_detrend_1D(self):
        x = np.arange(5)
        np.testing.assert_allclose(detrend_1D(1 + 2 * x), [0., 0., 0., 0., 0.])
        np.testing.assert_allclose(
            detrend_1D(x**0.5),
            [-0.28284271,  0.24395221,  0.18496069,  0.02959285, -0.17566304])

    def test_linregress(self):
        x = np.arange(3*4).reshape((3, 4))
        y = 3 + 2 * x
        np.testing.assert_allclose(
            linregress(x, y, axis=0),
            ([[3., 3., 3., 3.]],
             [[2., 2., 2., 2.]],
             [[1., 1., 1., 1.]]))
        np.testing.assert_allclose(
            linregress(x, y, axis=0),
            ([[3., 3., 3., 3.]],
             [[2., 2., 2., 2.]],
             [[1., 1., 1., 1.]]))
        np.testing.assert_allclose(
            linregress(x, y, axis=1),
            ([[3.], [3.], [3.]],
             [[2.], [2.], [2.]],
             [[1.], [1.], [1.]]))
        with self.assertRaises(np.AxisError):
            linregress(x, y, axis=2)

    def test_detrend(self):
        x = np.arange(3*4).reshape((3, 4))
        y = 3 + 2 * x
        np.testing.assert_allclose(detrend(y, axis=0), 0.)
        np.testing.assert_allclose(detrend(y, axis=1), 0.)
        with self.assertRaises(IndexError):
            detrend(y, axis=2)

unittest.main(argv=[''], verbosity=2, exit=False)

test_detrend (__main__.TestNotebook) ... ok
test_detrend_1D (__main__.TestNotebook) ... ok
test_linregress (__main__.TestNotebook) ... ok
test_linregress_1D (__main__.TestNotebook) ... ok

----------------------------------------------------------------------
Ran 4 tests in 0.004s

OK

<unittest.main.TestProgram at 0x7f7a26f52bf0>

Doctests

import doctest

doctest.testmod(verbose=True)

Trying:
    x = np.arange(3*4).reshape((3, 4))
Expecting nothing
ok
Trying:
    y = 3 + 2 * x
Expecting nothing
ok
Trying:
    detrend(y, axis=0)  # doctest: +NORMALIZE_WHITESPACE
Expecting:
    array([[0., 0., 0., 0.],
           [0., 0., 0., 0.],
           [0., 0., 0., 0.]])
ok
Trying:
    detrend(y, axis=1)  # doctest: +NORMALIZE_WHITESPACE
Expecting:
    array([[0., 0., 0., 0.],
           [0., 0., 0., 0.],
           [0., 0., 0., 0.]])
ok
Trying:
    detrend(y, axis=2)
Expecting:
    Traceback (most recent call last):
        ...
    IndexError: index 2 is out of bounds for axis 0 with size 2
ok
Trying:
    x = np.arange(5)
Expecting nothing
ok
Trying:
    detrend_1D(1 + 2 * x)
Expecting:
    array([0., 0., 0., 0., 0.])
ok
Trying:
    detrend_1D(x**0.5)
Expecting:
    array([-0.28284271,  0.24395221,  0.18496069,  0.02959285, -0.17566304])
ok
Trying:
    x = np.arange(3*4).reshape((3, 4))
Expecting nothing
ok
Trying:
    y = 3 + 2 * x
Expecting nothing
ok
Trying:
    linregress(x, y, axis=0)  # doctest: +NORMALIZE_WHITESPACE
Expecting:
    (array([[3., 3., 3., 3.]]),
     array([[2., 2., 2., 2.]]),
     array([[1., 1., 1., 1.]]))
ok
Trying:
    linregress(x, y, axis=1)  # doctest: +NORMALIZE_WHITESPACE
Expecting:
    (array([[3.], [3.], [3.]]),
     array([[2.], [2.], [2.]]),
     array([[1.], [1.], [1.]]))
ok
Trying:
    linregress(x, y, axis=2)
Expecting:
    Traceback (most recent call last):
        ...
    numpy.AxisError: axis 2 is out of bounds for array of dimension 2
ok
Trying:
    x = np.arange(10)
Expecting nothing
ok
Trying:
    linregress_1D(x, 3 + 2 * x)
Expecting:
    (3.0, 2.0, 1.0)
ok
Trying:
    linregress_1D(x, 3 + 2 * x + np.round(np.sin(x * 10), 3))
Expecting:
    (2.8803090909090923, 2.0317757575757573, 0.9925669351010524)
ok
6 items had no tests:
    __main__
    __main__.TestNotebook
    __main__.TestNotebook.test_detrend
    __main__.TestNotebook.test_detrend_1D
    __main__.TestNotebook.test_linregress
    __main__.TestNotebook.test_linregress_1D
4 items passed all tests:
   5 tests in __main__.detrend
   3 tests in __main__.detrend_1D
   5 tests in __main__.linregress
   3 tests in __main__.linregress_1D
16 tests in 10 items.
16 passed and 0 failed.
Test passed.

TestResults(failed=0, attempted=16)

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