Les valeurs observées obtenues par des expériences sont généralement des valeurs discrètes, mais il y a des moments où vous voulez trouver des valeurs entre les deux. À ce moment-là, diverses méthodes d'interpolation (pas de complémentation) et d'autres ajustements de courbe (approximation de courbe) sont utilisées. Ces solutions sont souvent utilisées comme sujet d'exercices de programmation, mais il existe de nombreuses bibliothèques puissantes dans la pratique, il est donc plus pratique d'utiliser une bibliothèque existante que de créer la vôtre.

Expérience informatique

Ici, les trois fonctions mathématiques suivantes f (x), g (x) et h (x) sont préparées sous forme d'expériences informatiques.

import numpy as np
f = lambda x: 1/(1 + np.exp(-x))
g = lambda x: 1.0/(1.0+x**2)
h = lambda x: np.sin(x)

On suppose que les trois fonctions mathématiques ci-dessus montrent le phénomène qui se produit réellement.

Observations expérimentales

Cependant, en réalité, les points d'observation obtenus dans l'expérience sont des valeurs discrètes comme suit.

x_observed = np.linspace(-10, 10, 11)
x_observed

array([-10.,  -8.,  -6.,  -4.,  -2.,   0.,   2.,   4.,   6.,   8.,  10.])

À ce stade, illustrons le type de valeur observée par l'expérience.

%matplotlib inline
import matplotlib.pyplot as plt

for func in [f, g, h]:
    y_observed = func(x_observed)
    plt.scatter(x_observed, y_observed)
    plt.grid()
    plt.show()

Quel genre de courbe pouvez-vous imaginer à partir des observations ci-dessus?

Vérité inobservée

Illustrons la vérité qui n'est pas observée dans l'expérience.

x_latent = np.linspace(-10, 10, 101)
x_latent

array([-10. ,  -9.8,  -9.6,  -9.4,  -9.2,  -9. ,  -8.8,  -8.6,  -8.4,
        -8.2,  -8. ,  -7.8,  -7.6,  -7.4,  -7.2,  -7. ,  -6.8,  -6.6,
        -6.4,  -6.2,  -6. ,  -5.8,  -5.6,  -5.4,  -5.2,  -5. ,  -4.8,
        -4.6,  -4.4,  -4.2,  -4. ,  -3.8,  -3.6,  -3.4,  -3.2,  -3. ,
        -2.8,  -2.6,  -2.4,  -2.2,  -2. ,  -1.8,  -1.6,  -1.4,  -1.2,
        -1. ,  -0.8,  -0.6,  -0.4,  -0.2,   0. ,   0.2,   0.4,   0.6,
         0.8,   1. ,   1.2,   1.4,   1.6,   1.8,   2. ,   2.2,   2.4,
         2.6,   2.8,   3. ,   3.2,   3.4,   3.6,   3.8,   4. ,   4.2,
         4.4,   4.6,   4.8,   5. ,   5.2,   5.4,   5.6,   5.8,   6. ,
         6.2,   6.4,   6.6,   6.8,   7. ,   7.2,   7.4,   7.6,   7.8,
         8. ,   8.2,   8.4,   8.6,   8.8,   9. ,   9.2,   9.4,   9.6,
         9.8,  10. ])

fig_idx = 0
plt.figure(figsize=(12,4))
for func in [f, g, h]:
    fig_idx += 1
    y_observed = func(x_observed)
    y_latent = func(x_latent)
    plt.subplot(1, 3, fig_idx)
    plt.scatter(x_observed, y_observed)
    plt.plot(x_latent, y_latent)
    plt.grid()
plt.show()

Maintenant, ce que je veux demander ici, c'est à quel point ce "vrai chiffre" peut être déduit des valeurs d'observation limitées.

Interpolation avec Scipy.interpolate

scipy fournit une bibliothèque puissante appelée interpolation qui vous permet d'utiliser une variété de méthodes d'interpolation.

from scipy import interpolate

ip1 = ["Interpolation du point le plus proche", lambda x, y: interpolate.interp1d(x, y, kind="nearest")]
ip2 = ["Interpolation linéaire", interpolate.interp1d]
ip3 = ["Interpolation de Lagrange", interpolate.lagrange]
ip4 = ["Interpolation du centre de gravité", interpolate.BarycentricInterpolator]
ip5 = ["Interpolation Krogh", interpolate.KroghInterpolator]
ip6 = ["Interpolation spline de second ordre", lambda x, y: interpolate.interp1d(x, y, kind="quadratic")]
ip7 = ["Interpolation spline d'ordre 3", lambda x, y: interpolate.interp1d(x, y, kind="cubic")]
ip8 = ["Interpolation d'automne", interpolate.Akima1DInterpolator]
ip9 = ["Interpolation Hermite cubique sectionnelle", interpolate.PchipInterpolator]

x_observed = np.linspace(-10, 10, 11)
x_observed

array([-10.,  -8.,  -6.,  -4.,  -2.,   0.,   2.,   4.,   6.,   8.,  10.])

y_observed = f(x_observed)
y_latent = f(x_latent)
for method_name, method in [ip1, ip2, ip3, ip4, ip5, ip6, ip7, ip8, ip9]:
    print(method_name)
    fitted_curve = method(x_observed, y_observed)
    plt.scatter(x_observed, y_observed, label="observed")
    plt.plot(x_latent, fitted_curve(x_latent), c="red", label="fitted")
    plt.plot(x_latent, y_latent, label="latent")
    plt.grid()
    plt.legend()
    plt.show()

Interpolation du point le plus proche

Interpolation linéaire

Interpolation de Lagrange

Interpolation du centre de gravité

Interpolation Krogh

Interpolation spline de second ordre

Interpolation spline d'ordre 3

Interpolation d'automne

Interpolation Hermite cubique sectionnelle

Défi 1

Interpoler g (x) avec diverses méthodes d'interpolation.
Interpoler h (x) avec diverses méthodes d'interpolation.
Interpole les mesures suivantes avec diverses méthodes d'interpolation.

x_observed = np.array([9, 28, 38, 58, 88, 98, 108, 118, 128, 138, 148, 158, 168, 178, 188, 198, 208, 218, 228, 238, 278, 288, 298])
y_observed = np.array([51, 80, 112, 294, 286, 110, 59, 70, 56, 70, 104, 59, 59, 72, 87, 99, 64, 60, 74, 151, 157, 57, 83])

x_latent = np.linspace(min(x_observed), max(x_observed), 100)

Ajustement de courbe avec Numpy.polyfit

.Polfyfit de Numpy facilite l'ajustement des courbes avec la méthode des moindres carrés. Si vous recalculez les valeurs d'observation sautantes observées par l'expérience

x_observed = np.linspace(-10, 10, 11)
x_observed

array([-10.,  -8.,  -6.,  -4.,  -2.,   0.,   2.,   4.,   6.,   8.,  10.])

y_observed = f(x_observed)
y_observed

array([4.53978687e-05, 3.35350130e-04, 2.47262316e-03, 1.79862100e-02,
       1.19202922e-01, 5.00000000e-01, 8.80797078e-01, 9.82013790e-01,
       9.97527377e-01, 9.99664650e-01, 9.99954602e-01])

Lorsque la régression linéaire est effectuée sur les x et y ci-dessus par la méthode du carré minimum,

coefficients = np.polyfit(x_observed, y_observed, 1)
coefficients

array([0.06668944, 0.5       ])

La valeur de retour ci-dessus correspond directement aux coefficients a et b de y = ax + b. Exprimé avec Sympy

import sympy as sym
from sympy.plotting import plot
sym.init_printing(use_unicode=True)
%matplotlib inline

x, y = sym.symbols("x y")

#Si vous utilisez Google Colab, exécutez pour prendre en charge l'affichage TeX par Sympy
def custom_latex_printer(exp,**options):
    from google.colab.output._publish import javascript
    url = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.3/latest.js?config=default"
    javascript(url=url)
    return sym.printing.latex(exp,**options)

sym.init_printing(use_latex="mathjax",latex_printer=custom_latex_printer)

expr = 0
for index, coefficient in enumerate(coefficients):
    expr += coefficient * x ** (len(coefficients) - index - 1)
display(sym.Eq(y, expr))

y = 0.0666894399883493 x + 0.5

C'est l'équation de régression obtenue.

Le dernier argument vous permet de spécifier l'ordre de retour. Par exemple, si vous souhaitez intégrer une formule cubique

coefficients = np.polyfit(x_observed, y_observed, 3)

expr = 0
for index, coefficient in enumerate(coefficients):
    expr += coefficient * x ** (len(coefficients) - index - 1)
display(sym.Eq(y, expr))

y = - 0.000746038674650085 x^{3} + 0.119807393623435 x + 0.5

Cela montre que nous avons pu revenir.

Pour dessiner la courbe résultante, nous avons besoin de plusieurs x et de leurs y correspondants, que nous obtenons:

fitted_curve = np.poly1d(np.polyfit(x_observed, y_observed, 3))(x_latent)
fitted_curve

array([ 0.04796474,  0.02805317,  0.00989629, -0.00654171, -0.02129666,
       -0.03440435, -0.0459006 , -0.05582121, -0.064202  , -0.07107878,
       -0.07648735, -0.08046353, -0.08304312, -0.08426194, -0.08415579,
       -0.08276049, -0.08011184, -0.07624566, -0.07119776, -0.06500394,
       -0.05770001, -0.04932179, -0.03990508, -0.02948569, -0.01809944,
       -0.00578213,  0.00743042,  0.02150241,  0.03639803,  0.05208146,
        0.0685169 ,  0.08566854,  0.10350056,  0.12197717,  0.14106254,
        0.16072086,  0.18091634,  0.20161315,  0.22277549,  0.24436755,
        0.26635352,  0.28869759,  0.31136394,  0.33431678,  0.35752028,
        0.38093865,  0.40453606,  0.42827671,  0.45212479,  0.47604449,
        0.5       ,  0.52395551,  0.54787521,  0.57172329,  0.59546394,
        0.61906135,  0.64247972,  0.66568322,  0.68863606,  0.71130241,
        0.73364648,  0.75563245,  0.77722451,  0.79838685,  0.81908366,
        0.83927914,  0.85893746,  0.87802283,  0.89649944,  0.91433146,
        0.9314831 ,  0.94791854,  0.96360197,  0.97849759,  0.99256958,
        1.00578213,  1.01809944,  1.02948569,  1.03990508,  1.04932179,
        1.05770001,  1.06500394,  1.07119776,  1.07624566,  1.08011184,
        1.08276049,  1.08415579,  1.08426194,  1.08304312,  1.08046353,
        1.07648735,  1.07107878,  1.064202  ,  1.05582121,  1.0459006 ,
        1.03440435,  1.02129666,  1.00654171,  0.99010371,  0.97194683,
        0.95203526])

plt.scatter(x_observed, f(x_observed), label="observed")
plt.plot(x_latent, fitted_curve, c="red", label="fitted")
plt.plot(x_latent, f(x_latent), label="latent")
plt.grid()
plt.legend()

<matplotlib.legend.Legend at 0x7fb673a3c630>

Défi 2

Approximative g (x) avec une équation quintique.
Approximation de h (x) avec une équation d'ordre 9.
Approximer les valeurs mesurées suivantes avec une équation d'ordre 9.

x_observed = np.array([9, 28, 38, 58, 88, 98, 108, 118, 128, 138, 148, 158, 168, 178, 188, 198, 208, 218, 228, 238, 278, 288, 298])
y_observed = np.array([51, 80, 112, 294, 286, 110, 59, 70, 56, 70, 104, 59, 59, 72, 87, 99, 64, 60, 74, 151, 157, 57, 83])

Ajustement de courbe avec Scipy.optimize.curve_fit

Une façon de faire un ajustement de courbe (approximation de courbe) à l'aide de Python est d'utiliser scipy.optimize.curve_fit.

Valeurs expérimentales utilisées

x_observed = [9, 28, 38, 58, 88, 98, 108, 118, 128, 138, 148, 158, 168, 178, 188, 198, 208, 218, 228, 238, 278, 288, 298]
y_observed = [51, 80, 112, 294, 286, 110, 59, 70, 56, 70, 104, 59, 59, 72, 87, 99, 64, 60, 74, 151, 157, 57, 83]

Passons du type List de Python au type Array de Numpy.

import numpy as np
x_observed = np.array(x_observed)
y_observed = np.array(y_observed)

Approximation linéaire

Commençons par l'approximer à une équation linéaire. Définissez la fonction func1. a et b sont les valeurs que vous recherchez.

def func1(X, a, b): #Approximation linéaire
    Y = a + b * X
    return Y

Voici un exemple d'utilisation de la fonction func1.

func1(x_observed, 2, 2)

array([ 20,  58,  78, 118, 178, 198, 218, 238, 258, 278, 298, 318, 338,
       358, 378, 398, 418, 438, 458, 478, 558, 578, 598])

Maintenant, faisons une approximation en utilisant scipy.optimize.curve_fit.

from scipy.optimize import curve_fit  
popt, pcov = curve_fit(func1,x_observed,y_observed) #popt est l'estimation optimale, pcov est la covariance
popt

array([125.77023172,  -0.1605313 ])

Le popt obtenu ici stocke l'estimation optimale. Comparons-le avec l'estimation obtenue par ajustement de courbe à l'aide de Numpy.polyfit.

Approximation quadratique

Ensuite, faisons une approximation à une équation quadratique. Définissez la fonction func2. a, b et c sont les valeurs souhaitées.

def func2(X, a, b, c): #Approximation quadratique
    Y = a + b * X + c * X ** 2
    return Y

Voici un exemple d'utilisation de la fonction func2.

func2(x_observed, 2, 2, 2)

array([   182,   1626,   2966,   6846,  15666,  19406,  23546,  28086,
        33026,  38366,  44106,  50246,  56786,  63726,  71066,  78806,
        86946,  95486, 104426, 113766, 155126, 166466, 178206])

Maintenant, faisons une approximation en utilisant scipy.optimize.curve_fit.

from scipy.optimize import curve_fit  
popt, pcov = curve_fit(func2,x_observed,y_observed) #popt est l'estimation optimale, pcov est la covariance
popt

array([ 1.39787684e+02, -4.07809516e-01,  7.97183226e-04])

Le popt obtenu ici stocke l'estimation optimale. Comparons-le avec l'estimation obtenue par ajustement de courbe à l'aide de Numpy.polyfit.

En procédant comme suit, vous pouvez dessiner une courbe approximative en utilisant la valeur estimée optimale.

func2(x_observed, 1.39787684e+02, -4.07809516e-01,  7.97183226e-04)

array([136.1819702 , 128.9940092 , 125.44205497, 118.81645644,
       110.07383349, 107.47849913, 105.04260142, 102.76614035,
       100.64911593,  98.69152815,  96.89337701,  95.25466253,
        93.77538468,  92.45554348,  91.29513893,  90.29417102,
        89.45263976,  88.77054514,  88.24788717,  87.88466585,
        88.02614699,  88.46010889,  89.05350743])

Approximation cubique

De même, approchons-le à une équation cubique. Définissez la fonction func3. a, b, c et d sont les valeurs souhaitées.

def func3(X, a, b, c, d): #Approximation cubique
    Y = a + b * X + c * X ** 2 + d * X ** 3
    return Y

Maintenant, faisons une approximation en utilisant scipy.optimize.curve_fit.

from scipy.optimize import curve_fit  
popt, pcov = curve_fit(func3,x_observed,y_observed) #popt est l'estimation optimale, pcov est la covariance
popt

array([ 7.84214107e+01,  1.88213104e+00, -1.74165777e-02,  3.89638087e-05])

Le popt obtenu ici stocke l'estimation optimale. Comparons-le avec l'estimation obtenue par ajustement de courbe à l'aide de Numpy.polyfit.

En procédant comme suit, vous pouvez dessiner une courbe approximative en utilisant la valeur estimée optimale.

func3(x_observed, 7.84214107e+01,  1.88213104e+00, -1.74165777e-02,  3.89638087e-05)

array([ 93.97825188, 118.32181643, 126.93087413, 136.59795028,
       135.72770915, 132.27386543, 127.62777811, 122.02323006,
       115.69400413, 108.87388316, 101.79665001,  94.69608754,
        87.80597859,  81.36010602,  75.59225267,  70.73620141,
        67.02573508,  64.69463654,  63.97668864,  65.10567422,
        92.76660851, 106.63700433, 123.75703075])

Je souhaite créer une fonction polyvalente utilisable quelle que soit la commande

Jusqu'à présent, j'ai créé des fonctions séparées pour les expressions linéaires, les expressions quadratiques et les expressions cubiques, mais c'est trop gênant, alors créons une fonction à usage général qui peut être utilisée quel que soit le nombre d'expressions. L'ordre d'approximation est automatiquement déterminé par le nombre de paramètres.

import numpy as np
def func(X, *params):
    Y = np.zeros_like(X)
    for i, param in enumerate(params):
        Y = Y + np.array(param * X ** i)
    return Y

Exécutez les deux codes suivants pour vérifier que le calcul func2 ci-dessus est le même que le calcul func3.

func(x_observed, 1.39787684e+02, -4.07809516e-01,  7.97183226e-04)

array([136.1819702 , 128.9940092 , 125.44205497, 118.81645644,
       110.07383349, 107.47849913, 105.04260142, 102.76614035,
       100.64911593,  98.69152815,  96.89337701,  95.25466253,
        93.77538468,  92.45554348,  91.29513893,  90.29417102,
        89.45263976,  88.77054514,  88.24788717,  87.88466585,
        88.02614699,  88.46010889,  89.05350743])

func(x_observed, 7.84214107e+01,  1.88213104e+00, -1.74165777e-02,  3.89638087e-05)

array([ 93.97825188, 118.32181643, 126.93087413, 136.59795028,
       135.72770915, 132.27386543, 127.62777811, 122.02323006,
       115.69400413, 108.87388316, 101.79665001,  94.69608754,
        87.80597859,  81.36010602,  75.59225267,  70.73620141,
        67.02573508,  64.69463654,  63.97668864,  65.10567422,
        92.76660851, 106.63700433, 123.75703075])

Approximation polygonale avec une fonction à usage général qui peut être utilisée avec n'importe quelle expression de degré

À ce stade, une approximation polynomiale peut être effectuée sur toute expression de degré. Cependant, vous devrez entrer le nombre requis de valeurs initiales avec p0 = pour spécifier le nombre de paramètres.

from scipy.optimize import curve_fit  
popt, pcov = curve_fit(func,x_observed,y_observed, p0=[1, 1]) 
popt

array([125.77023172,  -0.1605313 ])

from scipy.optimize import curve_fit  
popt, pcov = curve_fit(func,x_observed,y_observed, p0=[1, 1, 1]) 
popt

array([ 1.39787684e+02, -4.07809516e-01,  7.97183226e-04])

from scipy.optimize import curve_fit  
popt, pcov = curve_fit(func,x_observed,y_observed, p0=[1, 1, 1, 1]) 
popt

array([ 7.84214107e+01,  1.88213104e+00, -1.74165777e-02,  3.89638087e-05])

from scipy.optimize import curve_fit  
popt, pcov = curve_fit(func,x_observed,y_observed, p0=[1, 1, 1, 1, 1]) 
popt

array([-7.54495613e+01,  1.13179580e+01, -1.50591743e-01,  7.02970546e-04,
       -1.07313869e-06])

Comparaison avec l'ajustement de courbe avec Numpy.polyfit

Comparons le résultat du calcul ci-dessus avec la valeur estimée obtenue par ajustement de courbe à l'aide de Numpy.polyfit.

Numpy.polyfit est pratique et facile à utiliser si vous voulez juste faire une approximation polymorphe, mais vous ne pouvez faire qu'une approximation polymorphe. Si vous souhaitez effectuer une approximation avec d'autres fonctions, il est préférable de comprendre comment utiliser scipy.optimize.curve_fit.

Je veux savoir quelles autres méthodes d'optimisation sont disponibles

Si vous dites "optimiser?", Il existe diverses autres méthodes expliquées. Bien que ce soit en anglais.

from scipy import optimize
optimize?

Bien sûr, vous pouvez utiliser Google, mais si vous utilisez "?", Vous pouvez accéder directement à l'explication, et il est bon que vous puissiez l'utiliser hors ligne.

Méthode de Newton et dichotomie

La dichotomie et la méthode de Newton sont des méthodes typiques pour trouver la solution de l'équation f (x) = 0. À titre d'exemple, utilisez cette fonction.

import math
def f(x):
    return math.exp(x) - 3 * x

import math
f = lambda x: math.exp(x) - 3 * x

Une bibliothèque pour le calcul scientifique et technologique appelée scipy est pratique.

from scipy import optimize

Dichotomie

La dichotomie peut être utilisée lorsque la fonction y = f (x) est continue dans l'intervalle [a, b] et qu'il existe une solution. Par exemple, si vous savez qu'il existe une solution dans l'intervalle [0,1]

optimize.bisect(f, 0, 1)

0.619061286737633

Oh facile.

Méthode Newton

La méthode de Newton peut être utilisée lorsque la fonction y = f (x) est monotone et continue sans point d'inflexion, et la dérivée de f (x) peut être obtenue. Lorsque vous voulez trouver x lorsque f (x) = 0

optimize.newton(f, 0)

0.6190612867359452

Oh facile.

Je souhaite en savoir plus sur l'utilisation

Vous pouvez utiliser "?" Pour vérifier les explications des paramètres et les exemples d'utilisation. Bien que ce soit en anglais.

optimize.bisect?

optimize.newton?

Optimisations telles que l'interpolation et l'ajustement de courbe

Expérience informatique

Observations expérimentales

Vérité inobservée

Interpolation avec Scipy.interpolate

Défi 1

Ajustement de courbe avec Numpy.polyfit

Défi 2

Ajustement de courbe avec Scipy.optimize.curve_fit

Valeurs expérimentales utilisées

Approximation linéaire

Approximation quadratique

Approximation cubique

Je souhaite créer une fonction polyvalente utilisable quelle que soit la commande

Approximation polygonale avec une fonction à usage général qui peut être utilisée avec n'importe quelle expression de degré

Comparaison avec l'ajustement de courbe avec Numpy.polyfit

Je veux savoir quelles autres méthodes d'optimisation sont disponibles

Méthode de Newton et dichotomie

Dichotomie

Méthode Newton

Je souhaite en savoir plus sur l'utilisation