Notebook originel: exam23_corrige.ipynb

Examen Science des données: dérivation numérique et propagation des incertitudes

Date: 14 novembre 2023, 15h45-17h45

Durée: 2 h

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

%matplotlib notebook

L’objectif général de l’exercice est la propagation des incertitudes, en utilisant les méthodes de dérivation numérique.

• La première partie concerne les fonctions univariées (\(\mathbb{R}\to\mathbb{R}\));

• la seconde partie aborde le cas bivarié (\(\mathbb{R}^{2}\to\mathbb{R}\)), en incluant la gestion des covariances;

• la troisième partie concerne le packaging des fonctions développées dans le cadre des deux premières parties.

Consignes

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

Cas des fonctions univariées (\(\mathbb{R}\to\mathbb{R}\))

Dérivation numérique

La dérivation numérique consiste à estimer les dérivées d’une fonction à l’aide de l’approximation des différences finies.

Ainsi, la dérivée première d’une fonction univariée (\(\mathbb{R}\to\mathbb{R}\)) peut s’approximer par le quotient de différence symétrique:

\[\begin{split}\begin{align} f'(x) &= \frac{f(x + h) - f(x - h)}{2h} + \mathcal{O}(h^{2}) \\ &\approx \frac{1}{h}\left[ w_{-1} f(x_{-1}) + w_{+1} f(x_{+1}) \right] \end{align}\end{split}\]

avec \(w_{\pm 1} = \pm 1/2\) et \(x_{\pm 1} = x \pm h\). Puisque l’approximation est en \(h^{2}\), la méthode est dite d’ordre deux (toutes les méthodes symétriques sont par construction d’ordre paire).

1. Écrire une fonction numdiff_univariate_D1_O2(func, x, h=1e-3) qui retourne la dérivée première en x de la fonction univariée func en utilisant la méthode du quotient de différence symétrique d’ordre 2.

def numdiff_univariate_D1_O2(func, x, h=1e-3):

    return (func(x + h) - func(x - h)) / (2 * h)

numdiff_univariate_D1_O2(lambda x: x**3, np.pi), 3 * np.pi ** 2

(29.608814203264444, 29.608813203268074)

La méthode précédente se généralise à d’autres ordres de précision en utilisant les coefficients appropriés (voir Finite Difference Coefficients), p.ex.

\[\begin{split}\begin{align} f'(x) &= \frac{1}{h}\left[1/12 f(x_{-2}) - 2/3 f(x_{-1}) + 2/3 f(x_{+1}) - 1/12 f(x_{+2})\right] + \mathcal{O}(h^{4}) \\ &\approx \frac{1}{h} \sum_{i=-2}^{+2} w_{i} f(x_{i}) \end{align}\end{split}\]

avec \(x_{i} = x + i h\) et les coefficients \(w_{i}\) donnés dans numdiff_D1_coeffs[4].

# Coefficients pour les différences centrales (source: https://en.wikipedia.org/wiki/Finite_difference_coefficient)
# numdiff_D1_coeffs = {ordre: coeffs}
numdiff_D1_coeffs = {2: (-1/2, 0, 1/2),
                     4: (1/12, -2/3, 0, 2/3, -1/12)}

2. Écrire une fonction numdiff_univariate_D1_O4(func, x, h=1e-3) retournant la dérivée première en x de la fonction univariée func en utilisant la méthode du quotient de différence symétrique d’ordre 4.

def numdiff_univariate_D1_O4(func, x, h=1e-3):

    coeffs = numdiff_D1_coeffs[4]

    s = 0
    for i in range(-2, 2+1):
        if coeffs[i + 2]:  # no need to compute if coeff = 0
            s += coeffs[i + 2] * func(x + i * h)

    return s / h

numdiff_univariate_D1_O4(lambda x: x**3, np.pi / 2), 3 * (np.pi / 2)**2

(7.40220330081559, 7.4022033008170185)

3. Écrire une fonction numdiff_univariate_D1(func, x, order=2, h=1e-3) généralisant la fonction précédente pour calculer la dérivée première à l’ordre \(k\) (\(k=2, 4\)) en utilisant les coefficients de numdiff_D1_coeffs. La fonction devra lever l’exception NotImplementedError pour les ordres non supportés.

def numdiff_univariate_D1(func, x, order=2, h=1e-3):

    if order in numdiff_D1_coeffs:
        coeffs = numdiff_D1_coeffs[order]
    else:
        raise NotImplementedError(f"1st derivative: order {order} not supported.")

    horder = order // 2
    s = 0
    for i in range(-horder, horder+1):
        if coeffs[i + horder]:  # no need to compute if coeff = 0
            s += coeffs[i + horder] * func(x + i * h)

    return s / h

4. En utilisant la fonction précédente, tracez la dérivée première de la fonction sin aux ordres 2 et 4 pour \(h=10^{-3}\) et \(h=10^{-6}\), et comparer à la solution analytique.

   

   Note: puisque la fonction numdiff_univariate_D1 n’est pas vectorisée (elle ne peut prendre qu’un x scalaire en entrée), utilisez une liste en compréhension pour calculer les différentes valeurs de la dérivée.

   Comme vous pouvez le constater, choisir un pas \(h\) arbitrairement petit n’est pas nécessairement plus précis, du fait des erreurs d’arrondis…

x = np.linspace(0, 2*np.pi)
func = np.sin
dfunc = np.cos

fig, (ax, ax2) = plt.subplots(2, 1)
ax.plot(x, dfunc(x), c='k', label="Analytic")
ax2.plot(x, np.zeros_like(x), c='k')
for order in (2, 4):
    for h in (1e-3, 1e-6):
        y = np.array([ numdiff_univariate_D1(func, x_, order=order, h=h) for x_ in x ])
        ax.plot(x, y, label=f"order={order}, h={h}")
        ax2.plot(x, np.abs(y / dfunc(x) - 1))
ax.set(xlabel="x", ylabel=f"{func.__name__}'", title=func.__name__)
ax.legend(fontsize='small')
ax2.set(xlabel="x", ylabel=f"|{func.__name__}' / {dfunc.__name__} - 1|", yscale='log');

Propagation des incertitudes

Il est possible de calculer l’écart standard attendu d’une expression mathématique à l’aide de l’approximation linéaire de la thérie de la propagation des incertitudes. Pour une fonction \(f\) (\(\mathbb{R}\to\mathbb{R}\)), l’incertidude \(\sigma_{f}\) sur \(f(x)\) associée à l’incertitude \(\sigma_{x}\) sur \(x\) s’écrit au premier ordre:

\[\sigma_{f} = \left|\frac{\text{d}f}{\text{d}x}\right| \sigma_{x}.\]

5. Écrire une fonction properr_univariate(f, x, dx, order=2, h=1e-3) calculant l’erreur sur \(f(x \pm dx)\) dans l’approximation linéaire.

def properr_univariate(f, x, dx):

    dfdx = numdiff_univariate_D1(f, x)

    return abs(dfdx) * dx

Cas bivarié (\(\mathbb{R}^{2}\to\mathbb{R}\))

La procédure précédente peut se généraliser aux fonctions multivariées (on se limitera ici au cas \(\mathbb{R}^{2} \to \mathbb{R}\)):

\[\sigma_{f}^{2} = (\partial_{x}f)^{2}\sigma_{x}^{2} + (\partial_{y}f)^{2}\sigma_{y}^{2} + 2 \partial_{x}f\,\partial_{y}f\,\text{Cov}(x, y)\]

où \(\text{Cov}(x, y) \triangleq \left\langle(x - \hat{x})(y - \hat{y})\right\rangle\) désigne la covariance entre \(x\) et \(y\), et \(\text{Corr}(x, y) = \text{Cov}(x, y) / (\sigma_{x}\sigma_{y})\) le coefficient de corrélation entre \(x\) et \(y\).

6. Écrire une fonction properr_bivariate(f, x, y, dx, dy, corrxy=0) permettant de propager les erreurs (et la corrélation) sur \(x\) et \(y\).

   Dans cette fonction, utiliser les fonctions partielles pour calculer les dérivées partielles nécessaires, p.ex.

   fx = lambda _x: f(_x, y)             # fonction partielle fx(x) = f(x, y) à y fixé
   dfdx = numdiff_univariate_D1(fx, x)  # df/dx|(x, y)
   

def properr_bivariate(f, x, y, dx, dy, corrxy=0):

    fx = lambda _x: f(_x, y)             # fonction partielle fx(x) = f(x, y) à y fixé
    dfdx = numdiff_univariate_D1(fx, x)  # df/dx|(x, y)

    fy = lambda _y: f(x, _y)             # fonction partielle fy(y) = f(x, y) à x fixé
    dfdy = numdiff_univariate_D1(fy, y)  # df/dy|(x, y)

    df2 = (dfdx * dx) ** 2 + (dfdy * dy) ** 2 + 2 * dfdx * dfdy * corrxy * dx * dy

    return np.abs(df2) ** 0.5

g = lambda x, y: np.sin(x * y**2)
x, dx = np.pi, 0.01
y, dy = 1, 0.1
print(g(x, y), properr_bivariate(g, x, y, dx, dy))
print(g(x, x), properr_bivariate(g, x, x, dx, dx, corrxy=1))

1.2246467991473532e-16 0.6283939694826167
-0.39828817883405304 0.27157691943916984

7. Vérifier la pertinence de l’utilisation de la covariance en calculant l’erreur sur \(g(x, x) = x^{2} - x^{2}\), p.ex. pour \(x = 1 \pm 0.1\).

g = lambda x, y: x**2 - y**2
x, dx = 1, 0.1
print(g(x, x), properr_bivariate(g, x, x, dx, dx))
print(g(x, y), properr_bivariate(g, x, x, dx, dx, corrxy=1))

0 0.2828427124745957
0 3.725290298461914e-09

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 et les inclure dans un fichier numdiff.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 Exam23_Nom_Prenom constitué du seul module numdiff.py, avec l’arborescence suivante:

   Exam23_Nom_Prenom/
   ├── numdiff_nom_prenom
   │   ├── numdiff.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. 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 Exam23_NOM_Prenom/
   

   et envoyer le tarball.

Tests

import unittest

class TestNotebook(unittest.TestCase):

    def setUp(self):

        self.f = lambda x: -x**3
        self.dfdx = lambda x: -3 * x**2
        self.x, self.dx = np.pi, 1e-2
        self.df = 0.2960881320326807  # uncertainties on f(x ± dx)
        self.y, self.dy = 1, 0.1
        self.g = lambda x, y: -np.sin(x * y**2)
        self.dg = 0.6283981031508404  # uncertainties on g(x ± dx, y ± dy)
        self.dgc = 0.2715898998939025  # uncertainties on g(x ± dx, x ± dx) (full correlation)
        self.places = 5

    def test_01_D1_O2(self):

        self.assertAlmostEqual(numdiff_univariate_D1_O2(self.f, self.x), self.dfdx(self.x), places=self.places)

    def test_02_D1_O4(self):

        self.assertAlmostEqual(numdiff_univariate_D1_O4(self.f, self.x), self.dfdx(self.x), places=self.places)

    def test_03_D1_value(self):

        self.assertAlmostEqual(numdiff_univariate_D1(self.f, self.x, order=2), self.dfdx(self.x),
                               places=self.places)
        self.assertAlmostEqual(numdiff_univariate_D1(self.f, self.x, order=4), self.dfdx(self.x),
                               places=self.places)

    def test_04_D1_except(self):

        with self.assertRaises(NotImplementedError):
            numdiff_univariate_D1(self.f, self.x, order=6)       # Invalid order

    def test_05_properr_u(self):

        self.assertAlmostEqual(properr_univariate(self.f, self.x, self.dx), self.df, places=self.places)
        self.assertAlmostEqual(properr_univariate(self.f, -self.x, self.dx), self.df, places=self.places)

    def test_06_properr_b(self):

        self.assertAlmostEqual(properr_bivariate(self.g, self.x, self.y, self.dx, self.dy), self.dg,
                               places=4)  # Uncorrelated case
        self.assertAlmostEqual(properr_bivariate(self.g, self.x, self.x, self.dx, self.dx, corrxy=1), self.dgc,
                               places=4)  # Fully correlated case
        self.assertAlmostEqual(properr_bivariate(lambda x, y: x**2 - y**2, 1, 1, 0.1, 0.1), 0.2828427124745957,
                               places=self.places)
        self.assertAlmostEqual(properr_bivariate(lambda x, y: x**2 - y**2, 1, 1, 0.1, 0.1, corrxy=1), 0,
                               places=self.places)

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

test_01_D1_O2 (__main__.TestNotebook) ... ok
test_02_D1_O4 (__main__.TestNotebook) ... ok
test_03_D1_value (__main__.TestNotebook) ... ok
test_04_D1_except (__main__.TestNotebook) ... ok
test_05_properr_u (__main__.TestNotebook) ... ok
test_06_properr_b (__main__.TestNotebook) ... ok

----------------------------------------------------------------------
Ran 6 tests in 0.004s

OK

<unittest.main.TestProgram at 0x7f04d4b139a0>

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