Python. Ajustement de courbe (curve fitting).

In [6]:
from termcolor import cprint
from dateutil.tz import gettz
import datetime as dt
import locale

locale.setlocale(locale.LC_ALL, ('fr_FR', 'UTF-8'))

tzi = gettz('Europe/Brussels')
cprint('{} {}'.format('danielhagnoul', dt.datetime.now(tz=tzi).isoformat()), 'cyan')
danielhagnoul 2020-04-29T23:48:06.703933+02:00

Wikipédia nous donne la définition suivante : "L'ajustement de courbe est une technique d'analyse d'une courbe expérimentale, consistant à construire une courbe à partir de fonctions mathématiques et d'ajuster les paramètres de ces fonctions pour se rapprocher de la courbe mesurée — on parle donc aussi d'ajustement de paramètres. On utilise souvent le terme anglais curve fitting, profile fitting ou simplement fitting, pour désigner cette méthode ; on utilise souvent le franglais « fitter une courbe » pour dire « ajuster une courbe ». (---) il s'agit (---) de régression polynomiale lorsque l'on utilise un polynôme pour simuler le phénomène (les paramètres physiques pouvant être déduits des coefficients du polynôme)."

Dans certains cas, on a un modèle théorique permettant de prévoir la forme de la courbe ; la méthode d'ajustement permet de déterminer les paramètres de l'échantillon. Dans d'autres cas, on utilise une fonction empirique ; on s'intéresse alors en général à la surface, la largeur ou à la position du maximum de la fonction. Pour les formes de type « pic » (distributions unimodales), on utilise fréquemment des fonctions gaussiennes, lorentziennes ou bien des combinaisons (fonctions ou pseudo-fonctions de Voigt), ou encore des fonctions de Pearson. Les courbes présentant un amortissement comportent fréquemment une composante exponentielle négative (fonction en e-x) ; les courbes en S peuvent être modélisées par une fonction sigmoïde.

Lorsqu'on ne dispose pas d'un modèle théorique, que l'on ne connait pas la forme de la courbe, on utiliser les méthodes numpy polifit et poly1d pour tracer la courbe à partir d'une équation polynomiale.

Exemple d'équation polynomiale

In [7]:
import numpy as np
import matplotlib.pyplot as plt

plt.title("Équation polynomiale y(t)")
plt.xlabel("x_data")
plt.ylabel("y_data")

x_data = np.array([1, 2, 3, 4, 5])
y_data = np.array([2, 0, 3, 7, 13])
plt.scatter(x_data, y_data, s=20, c="red", label="data")

y_params = np.polyfit(x_data, y_data, 3)  # équation de degré 3

""" Nota bene : pour trouver la courbe il faut faire plusieurs essais.
Le dégré 1 convient pour des points alignés sur une ligne.
Le degré 2 et plus pour les courbes plus ou moins complexes.
N'utiliser pas un degré si élevé qu'il force la courbe a passé par tous les points,
cette courbe ne représente pas bien vos données.
"""
y = np.poly1d(y_params)

print(y)

"""
      3         2
-0.25 x + 3.464 x - 10.29 x + 9
"""

t = np.linspace(min(x_data), max(x_data), 100)
plt.plot(t, y(t), c="blue", label="y(t)")

plt.legend()
plt.grid()
#plt.show()
3         2
-0.25 x + 3.464 x - 10.29 x + 9

Pour ajuster un ensemble de points à une équation, nous utiliseront la méthode scipy.optimize.curve_fit. Les arguments de la méthode curve_fit sont :

  • f (xdata, a, b, ...) : est la fonction d'ajustement où xdata est les données de la variable indépendante et a, b, ... sont les paramètres d'ajustement, aussi nombreux soient-ils, répertoriés comme arguments séparés. Évidemment, f (xdata, a, b, ...) doit retourner la valeur y de la fonction d'ajustement.
  • xdata : est le tableau contenant les données x.
  • ydata : est le tableau contenant les données y.
  • p0 : est un tuple contenant les suppositions initiales des paramètres d'ajustement. Les suppositions pour les paramètres d'ajustement sont définies égales à 1 si elles ne sont pas spécifiées. C'est presque toujours une bonne idée de spécifier les suppositions initiales pour les paramètres d'ajustement.
  • sigma : est le tableau contenant les incertitudes dans les données.
  • (** kwargs) : sont des arguments de mots-clés qui peuvent être passés à la routine d'ajustement scipy.optimize.leastsq que curve_fit appelle. Ceux-ci sont généralement laissés sans précision.

Exemple utilisant la méthode curve-fit

In [8]:
from scipy.optimize import curve_fit

def logistique(x, K, a, r):
    """ Fonction logistique (Verhulst)
    Paramètres K, a, r.
    Voir : https://fr.wikipedia.org/wiki/Fonction_logistique_(Verhulst)
    """
    y = K / (1.0 + a * np.exp(r * -x))
    return y


x = np.linspace(-10, 15, 100)
y = logistique(x, 150, 10, 0.65)

plt.plot(x, y, color="yellow", label="y")

""" Ajouter du bruit (Gaussian distribution) 
pour simuler les données expérimentales.
"""
noise = 1.10 * np.random.normal(size=y.size)
y_noise = y + noise

plt.scatter(x, y_noise, s=5, c="red", label="y_noise")

""" Curve Fit """
init_vals = [10, 1, 0.1]  # for [K, a, r]
best_vals, covar = curve_fit(logistique, x, y_noise, p0=init_vals)
print('best_vals: {}'.format(best_vals))

y_fit = logistique(x, best_vals[0], best_vals[1], best_vals[2])
plt.plot(x, y_fit, label="y_fit")

plt.legend()
plt.grid()
#plt.show()
best_vals: [149.81392941   9.69889248   0.64103518]

Exemple avec la loi de puissance

In [9]:
import matplotlib as mpl


def power_law(x, a, b):
    """ Calcul la loi de puissance avec les constantes a et b """
    return a*np.power(x, b)


# Générer un jeu de données factice
x_dummy = np.linspace(start=1, stop=1000, num=100)
y_dummy = power_law(x_dummy, 1, 0.5)

# Ajouter du bruit à partir d'une distribution gaussienne
noise = 1.5*np.random.normal(size=y_dummy.size)
y_dummy = y_dummy + noise

# Créer un objet figure et le stocker dans une variable appelée 'fig'
fig = plt.figure(figsize=(8, 5))

# Ajouter un objet axes à notre figure qui occupe la figure entière
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8])

# Modifier les graduations majeures et mineures des axes x et y
ax.xaxis.set_tick_params(which='major', size=10, width=2, direction='in')
ax.xaxis.set_tick_params(which='minor', size=7, width=2, direction='in')
ax.yaxis.set_tick_params(which='major', size=10,
                         width=2, direction='in', right='on')
ax.yaxis.set_tick_params(which='minor', size=7, width=2,
                         direction='in', right='on')


# Tracer les données exponentielles contenant du bruit (noise)
ax.scatter(x_dummy, y_dummy, s=20, color='#00b3b3', label='Data')

# Ajouter les étiquettes des axes x et y
ax.set_xlabel('x-data', labelpad=5)
ax.set_ylabel('y-data', labelpad=5)

# Définir la mise à l'échelle des axes x et y sur logarithmique
ax.set_xscale('log')
ax.set_yscale('log')

# Modifier les emplacements des graduations majeures et mineures des axes x et y
ax.xaxis.set_major_locator(mpl.ticker.LogLocator(base=10.0))
ax.yaxis.set_major_locator(mpl.ticker.LogLocator(base=10.0))

# Définissez les limites de l'axe
ax.set_xlim(10, 1000)
ax.set_ylim(1, 100)

# Ajuster les données factices de la loi de puissance
pars, cov = curve_fit(f=power_law, xdata=x_dummy, ydata=y_dummy, p0=[
                      0, 0], bounds=(-np.inf, np.inf))

print(pars)

# Tracer les données d'ajustement
ax.plot(x_dummy, power_law(x_dummy, *pars),
        linestyle='--', linewidth=2, color='black')

plt.grid()
#plt.show()
[1.05927151 0.49171581]
In [10]:
!jupyter nbconvert --to html curve_fitting.ipynb
[NbConvertApp] Converting notebook curve_fitting.ipynb to html
[NbConvertApp] Writing 302271 bytes to curve_fitting.html
In [ ]: