Beobachtete Werte, die durch Experimente erhalten wurden, sind normalerweise diskrete Werte, aber es gibt Zeiten, in denen Sie Werte dazwischen finden möchten. Zu diesem Zeitpunkt werden verschiedene Interpolationsmethoden (keine Komplementation) und andere Kurvenanpassungen (Kurvennäherung) verwendet. Diese Lösungen werden häufig als Gegenstand von Programmierübungen verwendet. In der Praxis gibt es jedoch viele leistungsstarke Bibliotheken. Daher ist es bequemer, eine vorhandene Bibliothek zu verwenden, als eine eigene zu erstellen.

Computerexperiment

Hier werden die folgenden drei mathematischen Funktionen f (x), g (x) und h (x) als Computerexperimente hergestellt.

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)

Es wird angenommen, dass die obigen drei mathematischen Funktionen das tatsächlich auftretende Phänomen zeigen.

Experimentelle Beobachtungen

In der Realität sind die im Experiment erhaltenen Beobachtungspunkte jedoch diskrete Werte wie folgt.

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

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

Lassen Sie uns zu diesem Zeitpunkt veranschaulichen, welche Art von Wert durch das Experiment beobachtet wird.

%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()

Welche Art von Kurve können Sie sich aus den obigen Beobachtungen vorstellen?

Unbeobachtete Wahrheit

Lassen Sie uns die Wahrheit veranschaulichen, die im Experiment nicht beobachtet wird.

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()

Nun möchte ich fragen, wie nahe diese "wahre Zahl" aus den begrenzten Beobachtungswerten abgeleitet werden kann.

Interpolation mit Scipy.interpolate

scipy bietet eine leistungsstarke Bibliothek namens Interpolate, mit der Sie verschiedene Interpolationsmethoden verwenden können.

from scipy import interpolate

ip1 = ["Nächste Punktinterpolation", lambda x, y: interpolate.interp1d(x, y, kind="nearest")]
ip2 = ["Lineare Interpolation", interpolate.interp1d]
ip3 = ["Lagrange-Interpolation", interpolate.lagrange]
ip4 = ["Schwerpunktinterpolation", interpolate.BarycentricInterpolator]
ip5 = ["Krogh-Interpolation", interpolate.KroghInterpolator]
ip6 = ["Spline-Interpolation zweiter Ordnung", lambda x, y: interpolate.interp1d(x, y, kind="quadratic")]
ip7 = ["Spline-Interpolation 3. Ordnung", lambda x, y: interpolate.interp1d(x, y, kind="cubic")]
ip8 = ["Herbstinterpolation", interpolate.Akima1DInterpolator]
ip9 = ["Schnittweise kubische Einsiedlerinterpolation", 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()

Nächste Punktinterpolation

Lineare Interpolation

Lagrange-Interpolation

Schwerpunktinterpolation

Krogh-Interpolation

Spline-Interpolation zweiter Ordnung

Spline-Interpolation 3. Ordnung

Herbstinterpolation

Schnittweise kubische Einsiedlerinterpolation

Herausforderung 1

Interpoliere g (x) mit verschiedenen Interpolationsmethoden.
Interpoliere h (x) mit verschiedenen Interpolationsmethoden.
Interpolieren Sie die folgenden Messungen mit verschiedenen Interpolationsmethoden.

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)

Kurvenanpassung mit Numpy.polyfit

Numpys .polfyfit macht es einfach, Kurvenanpassungen mit der Methode der kleinsten Quadrate durchzuführen. Wenn Sie die vom Experiment beobachteten Sprungbeobachtungswerte neu berechnen

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])

Wenn eine lineare Regression für die obigen x und y nach der Minimum-Square-Methode durchgeführt wird,

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

array([0.06668944, 0.5       ])

Der obige Rückgabewert entspricht direkt den Koeffizienten a und b von y = ax + b. Mit Sympy ausgedrückt

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

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

#Wenn Sie Google Colab verwenden, führen Sie es aus, um die TeX-Anzeige von Sympy zu unterstützen
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

Das ist die erhaltene Regressionsgleichung.

Mit dem letzten Argument können Sie die Reihenfolge der Rückgabe angeben. Zum Beispiel, wenn Sie in eine kubische Formel passen möchten

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

Es zeigt, dass wir zurückkehren konnten.

Um die resultierende Kurve zu zeichnen, benötigen wir mehrere x und die entsprechenden y, die wir erhalten:

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>

Herausforderung 2

Ungefähre g (x) mit einer Quintgleichung.
Ungefähre h (x) mit einer Gleichung 9. Ordnung.
Annähern Sie die folgenden Messwerte mit der Gleichung 9. Ordnung.

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])

Kurvenanpassung mit Scipy.optimize.curve_fit

Eine Möglichkeit zur Kurvenanpassung (Kurvenannäherung) mit Python ist die Verwendung von scipy.optimize.curve_fit.

Verwendete experimentelle Werte

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]

Lassen Sie uns vom Listentyp von Python in den Array-Typ von Numpy konvertieren.

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

Lineare Näherung

Lassen Sie es uns zunächst einer linearen Gleichung annähern. Definieren Sie die Funktion func1. a und b sind die Werte, die Sie finden möchten.

def func1(X, a, b): #Lineare Näherung
    Y = a + b * X
    return Y

Hier ist ein Beispiel für die Verwendung der Funktion 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])

Lassen Sie uns nun mit scipy.optimize.curve_fit approximieren.

from scipy.optimize import curve_fit  
popt, pcov = curve_fit(func1,x_observed,y_observed) #popt ist die optimale Schätzung, pcov ist die Kovarianz
popt

array([125.77023172,  -0.1605313 ])

Der hier erhaltene Popt speichert die optimale Schätzung. Vergleichen wir es mit der Schätzung, die durch Kurvenanpassung mit Numpy.polyfit erhalten wurde.

Quadratische Approximation

Als nächstes approximieren wir es einer quadratischen Gleichung. Definieren Sie die Funktion func2. a, b und c sind die gewünschten Werte.

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

Hier ist ein Beispiel für die Verwendung der Funktion 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])

Lassen Sie uns nun mit scipy.optimize.curve_fit approximieren.

from scipy.optimize import curve_fit  
popt, pcov = curve_fit(func2,x_observed,y_observed) #popt ist die optimale Schätzung, pcov ist die Kovarianz
popt

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

Der hier erhaltene Popt speichert die optimale Schätzung. Vergleichen wir es mit der Schätzung, die durch Kurvenanpassung mit Numpy.polyfit erhalten wurde.

Mit den folgenden Schritten können Sie eine ungefähre Kurve mit dem optimalen Schätzwert zeichnen.

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])

Kubische Annäherung

In ähnlicher Weise approximieren wir es einer kubischen Gleichung. Definieren Sie die Funktion func3. a, b, c und d sind die gewünschten Werte.

def func3(X, a, b, c, d): #Kubische Annäherung
    Y = a + b * X + c * X ** 2 + d * X ** 3
    return Y

Lassen Sie uns nun mit scipy.optimize.curve_fit approximieren.

from scipy.optimize import curve_fit  
popt, pcov = curve_fit(func3,x_observed,y_observed) #popt ist die optimale Schätzung, pcov ist die Kovarianz
popt

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

Der hier erhaltene Popt speichert die optimale Schätzung. Vergleichen wir es mit der Schätzung, die durch Kurvenanpassung mit Numpy.polyfit erhalten wurde.

Mit den folgenden Schritten können Sie eine ungefähre Kurve mit dem optimalen Schätzwert zeichnen.

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])

Ich möchte eine Allzweckfunktion erstellen, die unabhängig von der Reihenfolge verwendet werden kann

Bisher habe ich separate Funktionen für lineare Ausdrücke, quadratische Ausdrücke und kubische Ausdrücke erstellt. Dies ist jedoch zu problematisch. Erstellen Sie daher eine Allzweckfunktion, die unabhängig von der Anzahl der Ausdrücke verwendet werden kann. Die Reihenfolge der Approximation wird automatisch durch die Anzahl der Parameter bestimmt.

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

Führen Sie die folgenden zwei Codes aus, um sicherzustellen, dass die obige func2-Berechnung mit der func3-Berechnung übereinstimmt.

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])

Polygonale Approximation durch eine Allzweckfunktion, die unabhängig von der Reihenfolge verwendet werden kann

An diesem Punkt können Polynomnäherungen zu jedem Grad Ausdruck gemacht werden. Sie müssen jedoch die erforderliche Anzahl von Anfangswerten mit p0 = eingeben, um die Anzahl der Parameter anzugeben.

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])

Vergleich mit der Kurvenanpassung mit Numpy.polyfit

Vergleichen wir das obige Berechnungsergebnis mit dem geschätzten Wert, der durch Kurvenanpassung mit Numpy.polyfit erhalten wurde.

Numpy.polyfit ist bequem und einfach zu verwenden, wenn Sie nur eine polymorphe Approximation durchführen möchten, aber nur eine polymorphe Approximation durchführen können. Wenn Sie sich anderen Funktionen annähern möchten, ist es besser zu verstehen, wie scipy.optimize.curve_fit verwendet wird.

Ich möchte wissen, welche anderen Optimierungsmethoden verfügbar sind

Wenn Sie "Optimieren?" Sagen, werden verschiedene andere Methoden erläutert. Obwohl es auf Englisch ist.

from scipy import optimize
optimize?

Natürlich ist es in Ordnung zu googeln, aber wenn Sie "?" Verwenden, können Sie direkt zur Erklärung gehen, und es ist gut, dass Sie sie offline verwenden können.

Newton-Methode und Dichotomie

Die Dichotomie und die Newtonsche Methode sind typische Methoden, um die Lösung der Gleichung f (x) = 0 zu finden. Verwenden Sie als Beispiel diese Funktion.

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

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

Eine Bibliothek für wissenschaftliche und technologische Berechnungen namens scipy ist praktisch.

from scipy import optimize

Dichotomie

Die Dichotomie kann verwendet werden, wenn die Funktion y = f (x) im Intervall [a, b] stetig ist und es eine Lösung gibt. Zum Beispiel, wenn Sie wissen, dass es im Intervall [0,1] eine Lösung gibt.

optimize.bisect(f, 0, 1)

0.619061286737633

Oh einfach.

Newton-Methode

Die Newton-Methode kann verwendet werden, wenn die Funktion y = f (x) monoton und stetig ohne Wendepunkte ist und die Ableitung von f (x) erhalten werden kann. Wenn Sie x finden möchten, wenn f (x) = 0 ist

optimize.newton(f, 0)

0.6190612867359452

Oh einfach.

Ich möchte eine detailliertere Verwendung erfahren

Mit "?" Können Sie Parametererklärungen und Verwendungsbeispiele überprüfen. Obwohl es auf Englisch ist.

optimize.bisect?

optimize.newton?

[PYTHON] Optimierungen wie Interpolation und Kurvenanpassung

Computerexperiment

Experimentelle Beobachtungen

Unbeobachtete Wahrheit

Interpolation mit Scipy.interpolate

Herausforderung 1

Kurvenanpassung mit Numpy.polyfit

Herausforderung 2

Kurvenanpassung mit Scipy.optimize.curve_fit

Verwendete experimentelle Werte

Lineare Näherung

Quadratische Approximation

Kubische Annäherung

Ich möchte eine Allzweckfunktion erstellen, die unabhängig von der Reihenfolge verwendet werden kann

Polygonale Approximation durch eine Allzweckfunktion, die unabhängig von der Reihenfolge verwendet werden kann

Vergleich mit der Kurvenanpassung mit Numpy.polyfit

Ich möchte wissen, welche anderen Optimierungsmethoden verfügbar sind

Newton-Methode und Dichotomie

Dichotomie

Newton-Methode

Ich möchte eine detailliertere Verwendung erfahren