Approximieren Sie eine Bezier-Kurve durch einen bestimmten Punkt mit der Methode der kleinsten Quadrate in Python
Der Auslöser war, dass bei der Organisation der Programme, die ich in der Vergangenheit geschrieben habe, die Kurvenanpassung der Bezier-Kurve (Kurvenanpassung % 9A% E3% 81% 82% E3% 81% A6% E3% 81% AF% E3% 82% 81)) kam heraus, und ich fand es nostalgisch, aber die Formel war schwierig. Ich konnte es überhaupt nicht verstehen.
Warum kann ich es mit einer so mysteriösen Formel approximieren? Ich bin in den alten Tagen erstaunlich.
Anscheinend wurde es mit der Methode der kleinsten Quadrate angenähert, ich möchte es nach der Überprüfung in Python erklären. Ich zog das Nachschlagewerk heraus und studierte es erneut. Die Nachschlagewerke stammen übrigens von Dr. Kenichi Kanaya "Optimierte Mathematik, die Sie verstehen können - von Grundprinzipien bis zu Berechnungsmethoden". % E3% 81% 93% E3% 82% 8C% E3% 81% AA% E3% 82% 89% E5% 88% 86% E3% 81% 8B% E3% 82% 8B% E5% BF% 9C% E7 % 94% A8% E6% 95% B0% E5% AD% A6% E6% 95% 99% E5% AE% A4% E2% 80% 95% E6% 9C% 80% E5% B0% 8F% E4% BA % 8C% E4% B9% 97% E6% B3% 95% E3% 81% 8B% E3% 82% 89% E3% 82% A6% E3% 82% A7% E3% 83% BC% E3% 83% 96 % E3% 83% AC% E3% 83% 83% E3% 83% 88% E3% 81% BE% E3% 81% A7-% E9% 87% 91% E8% B0% B7-% E5% 81% A5 % E4% B8% 80 / dp / 4320017382). Ich empfehle dieses Buch, weil es sehr leicht zu verstehen ist.
Was ist die Methode der kleinsten Quadrate?
Minimum Square-Methode Ist eine typische Methode zur Approximation komplexer Daten und Funktionen und die wichtigste Grundlage für die Datenanalyse. Der Anwendungsbereich ist sehr breit und es handelt sich um eine Methode, die sowohl praktisch als auch schön ist und bei der jede Funktion angenähert werden kann, wenn sie differenziert werden kann.
Angenommen, Sie haben hier 20 Punkte.
Lassen Sie uns eine gerade Linie setzen, die den 20 Punkten am nächsten kommt.
Suchen Sie für die lineare Funktion $ f (x) = ax + b $ die Parameter $ a $, $ b $, die allen Punkten am nächsten kommen. Die typische Methode zum Finden von $ a $ und $ b $ ist die ** Methode der kleinsten Quadrate **.
Es wird als Methode der kleinsten Quadrate bezeichnet, da es die angegebene Position und die geplottete Position approximiert, so dass die Summe der Quadrate der ** Fehler ** minimiert wird **.
Gerade Annäherung
Versuchen wir zunächst die einfachste Näherung der geraden Linie und zeichnen eine gerade Linie, die sich den folgenden vier Punkten annähert.
points = [(107, 101), (449, 617), (816, 876), (1105,1153)]
Die lineare Funktion ist $ f (x) = ax + b $, $ a $ ist die Steigung und $ b $ ist der Abschnitt. Suchen Sie für die Funktion $ f (x) $ die Parameter $ a $ und $ b $, die das Quadrat des Fehlers am angegebenen Punkt minimieren.
Unter der Annahme, dass der angegebene Ort $ y $ ist, lautet die Funktion $ j (x, y) $, die das Quadrat des Fehlers berechnet, wie folgt.
j(x, y) = (y - (ax + b))^2
Da die Summe von $ j $ für N Punkte minimiert werden sollte, lautet die Formel für die Methode der kleinsten Quadrate wie folgt.
J = \frac {1}{2}\sum_{k=1}^{n} j(x_k, y_k)
Finden Sie die partielle Ableitung der Fehlerfunktion j (x, y)
Plötzlich kam das esoterische Wort ** Bias-Funktion ** heraus, aber mach dir um nichts Sorgen. Auch Grundschüler können rechnen, wenn sie die Berechnungsmethode kennen.
Die Differenzierung einer Funktion mit zwei oder mehr Variablen in Bezug auf eine Variable wird als ** partielle Differenzierung ** bezeichnet. Und die durch partielle Differenzierung erhaltene Funktion heißt ** partielle Ableitung **.
</ i> Erklärung für diejenigen, die die Unterscheidung vergessen haben
Die Differenzierung ist $ (X ^ n) '= nX ^ {n-1} $, und der Exponent ist n als Koeffizient und der Exponent ist n-1.
Beispiel:
Die Differenzierung ist die Steigung der Tangentenlinie des Diagramms, die so umformuliert werden kann, dass die momentane Geschwindigkeit des Diagramms ermittelt wird. Ein Punkt mit einer Differenz von 0 ist ein Wendepunkt im Fall eines konvexen Graphen.
Finden Sie das partielle Differential mit Sympy
Wenn Sie die Differenzierung verstehen, ist die teilweise Differenzierung nichts Neues, aber es ist in Ordnung, wenn Sie die schwierigen Dinge nicht verstehen. Wir haben Sympy (und Jupyter)!
Wenn Sie Sympy nicht kennen, lesen Sie bitte Erstellen Sie die stärkste Taschenrechnerumgebung mit Sympy + Jupyter.
Beginnen wir gleich mit Jupyter.
from sympy import *
init_session()
Damit ist das Sympy- und Jupyter-Setup abgeschlossen. Definieren Sie eine Funktion $ j $, um den Fehler zu finden.
a, b = symbols("a b")
j = (y - (a * x + b)) ** 2
Sie finden die partielle Differenzierung mit der Funktion sympy.diff. Zunächst wird die partielle Differenzierung für $ a $ durchgeführt.
Die Bias-Funktion für $ a $ wird als $ \ frac {\ partielle j} {\ partielle a} $ und für $ b $ als $ \ frac {\ partielle j} {\ partielle b} $ geschrieben.
j_a = diff(j, a)
j_b = diff(j, b)
\frac{\partial j}{\partial a}
=
- 2 x \left(- a x - b + y\right) \\
\frac{\partial j}{\partial b}
=
2 a x + 2 b - 2 y
Ich konnte die partielle Ableitung der linearen Funktion finden und mit "Sympy" wäre es ein einfacher Gewinn.
Die Methode des minimalen Quadrats beträgt 1/2 des Gesamtfehlers. Berechnen wir sie also.
Das Zuweisen eines Werts zu einem Ausdruck erfolgt mit der Funktion sympy.subs.
sum_a = sum([j_a.subs([(x, _x), (y, _y)]) for _x, _y in points]) / 2.
sum_b = sum([j_b.subs([(x, _x), (y, _y)]) for _x, _y in points]) / 2.
\frac {1}{2}\sum_{k=1}^{n} \frac{\partial j}{\partial a}
=
2099931.0 a + 2477.0 b - 2276721.0 \\
\frac {1}{2}\sum_{k=1}^{n} \frac{\partial j}{\partial b}
=
2477.0 a + 4 b - 2747.0
Sie haben die Hälfte jeder Summe und müssen nur die simultanen Gleichungen dieser beiden Gleichungen lösen.
\begin{bmatrix}
2099931.0 & 2477.0 \\
2477.0 & 4
\end{bmatrix}
\begin{bmatrix}
a \\
b
\end{bmatrix}
=
\begin{bmatrix}
2276721.0 \\
2747.0
\end{bmatrix}
solve([sum_a, sum_b], [a, b])
# {a:1.01694642025091,b:57.0059292596265}
Die Antwort lautet "a = 1.01694642025091", "b = 57.0059292596265". Setzen Sie es in die lineare Gleichung $ y = ax + b $ ein. Lassen Sie es uns gleich planen.
Da es sich um eine gerade Linie handelt, kann sie nicht alle Punkte durchlaufen, sondern einen nahe gelegenen Ort.
Bezierkurvenapproximation
Ich habe früher eine einfache gerade Linie angenähert, aber die Bezier-Kurve der kubischen Kurve kann durch genau das gleiche Verfahren angenähert werden.
In der Bezier-Kurve ist $ p_1 $ der Startpunkt, $ p_4 $ der Endpunkt und $ p_2 $ und $ p_3 $ die Kontrollpunkte.
Die Formel für die kubische Bezier-Kurve lautet:
bz = p_{1} \left(- t + 1\right)^{3} + p_{2} t \left(- 3 t + 3\right) \left(- t + 1\right) + p_{3} t^{2} \left(- 3 t + 3\right) + p_{4} t^{3}
Der Bereich von t ist 0 ≤ t ≤ 1.
Unter der Annahme, dass die Formel der Bezier-Kurve $ bz $ lautet und der zu approximierende Punkt $ p_x $ ist, lautet die Funktion $ j $ zur Berechnung des Fehlers wie folgt.
j = (p_x - bz)^{2}
Die Definition der Methode der kleinsten Quadrate war wie folgt.
J = \frac {1}{2}\sum_{k=1}^{n} j(x_k, y_k)
Die Parameter, die Sie suchen möchten, sind $ p_2 $ und $ p_3 $. Nehmen Sie daher mit jedem Parameter eine teilweise Differenzierung vor. Es gibt 1/2 der Definition der Methode der kleinsten Quadrate, also multiplizieren wir die Abweichungsfunktion im Voraus mit 1/2.
Schreiben wir sofort mit "Sympy".
p0, p1, p2, p3, p4, px, t = symbols("p0 p1 p2 p3 p4 px t")
bz = (1-t)**3*p1 + 3*(1-t)**2*t*p2 + 3*(1-t)*t**2*p3 + t**3*p4
j = (px - bz) ** 2
jp2 = simplify(diff(j, p2) / 2)
jp3 = simplify(diff(j, p3) / 2)
Die Hälfte der Bias-Funktionen von $ p_2 $ und $ p_3 $ sind wie folgt.
\frac{1}{2} \frac{\partial j}{\partial p_2}
=
3 t \left(t - 1\right)^{2} \left(- p_{1} \left(t - 1\right)^{3} + 3 p_{2} t \left(t - 1\right)^{2} - 3 p_{3} t^{2} \left(t - 1\right) + p_{4} t^{3} - px\right)
\\
\frac{1}{2} \frac{\partial j}{\partial p_3}
=
3 t^{2} \left(t - 1\right) \left(p_{1} \left(t - 1\right)^{3} - 3 p_{2} t \left(t - 1\right)^{2} + 3 p_{3} t^{2} \left(t - 1\right) - p_{4} t^{3} + px\right)
Definieren Sie den Paarwert des Punktes $ p_x $, den Sie übergeben möchten, und $ t $ zu diesem Zeitpunkt. Angenommen, Sie haben die folgenden 4 Punkte.
points = [(0, 0), (0.2, 65), (0.7, 45), (1.0, 100)]
Die Start- und Endpunkte sind $ p_1 $ und $ p_4 $, daher sind diese beiden Werte bekannt. Lösen Sie für $ p_2 $ und $ p_3 $, die Sie finden möchten, die Summe der Abweichungsfunktionen mit einer Reihengleichung.
const = ((p1, points[0][1]), (p4, points[-1][1]))
tp2 = sum([jp2.subs(const + ((t, x[0]), (px, x[1]))) for x in points[1:-1]])
tp3 = sum([jp3.subs(const + ((t, x[0]), (px, x[1]))) for x in points[1:-1]])
solve([tp2, tp3], [p2, p3])
# {p2:180.456349206349,p3:−53.0753968253968}
Die Antwort lautet "p2 = 180,456349206349, p3 = -53,0753968253968".
Sie können es richtig approximieren!
Schauen wir uns verschiedene Parameter an.
points = [(0, 0), (0.3, 50), (0.7, 80), (1.0, 100)]
points = [(0, 0), (0.2, 34), (0.4, 44), (0.6, 46), (0.8, 60), (1.0, 100)]
Sofern es sich nicht um einen sehr extremen Wert handelt, scheint er ordnungsgemäß zu bestehen. Natürlich gibt es einige Werte, die nicht übergeben werden können, aber irgendwie zieht es eine solche Linie.
Versuchen Sie, verschiedene Kurven zu erzeugen
Schließlich wird auf Jupyter die Bezier-Kurve durch die Methode des minimalen Quadrats angenähert, und die zu zeichnenden Zellen werden beschrieben.
%matplotlib inline
import matplotlib.pyplot as plt
p0, p1, p2, p3, p4, px, t = symbols("p0 p1 p2 p3 p4 px t")
bz = (1-t)**3*p1 + 3*(1-t)**2*t*p2 + 3*(1-t)*t**2*p3 + t**3*p4
j = (px - bz) ** 2
jp2 = simplify(diff(j, p2) / 2)
jp3 = simplify(diff(j, p3) / 2)
def plot(points):
const = ((p1, points[0][1]), (p4, points[-1][1]))
tp2 = sum([jp2.subs(const + ((t, x[0]), (px, x[1]))) for x in points[1:-1]])
tp3 = sum([jp3.subs(const + ((t, x[0]), (px, x[1]))) for x in points[1:-1]])
vec = solve([tp2, tp3], [p2, p3])
const = {p1: points[0][1], p2: vec[p2], p3: vec[p3], p4: points[-1][1]}
bz2 = bz.subs(const)
x = np.arange(0, 1+0.01, 0.01)
y = [
bz2.subs(t, _t)
for _t in x
]
plt.plot(x, y)
plt.scatter([p[0] for p in points], [p[1] for p in points])
Wenn Sie eine Liste von Punkten eingeben, die die Funktion "Plot" durchlaufen soll, bemüht sie sich, diese Punkte mit der Methode der kleinsten Quadrate zu durchlaufen. Die Verwendung wird wie folgt aufgerufen.
points = [(0, 0), (0.1, -40), (0.5, -50), (0.8, 30), (1.0, 100)]
plot(points)
Lösen Sie ohne Sympy
Schreiben Sie eine Implementierung in einer Umgebung ohne Sympy wie C / C ++.
\frac{1}{2} \frac{\partial j}{\partial p_2}
=
3 t \left(t - 1\right)^{2} \left(- p_{1} \left(t - 1\right)^{3} + 3 p_{2} t \left(t - 1\right)^{2} - 3 p_{3} t^{2} \left(t - 1\right) + p_{4} t^{3} - px\right)
\\
\frac{1}{2} \frac{\partial j}{\partial p_3}
=
3 t^{2} \left(t - 1\right) \left(p_{1} \left(t - 1\right)^{3} - 3 p_{2} t \left(t - 1\right)^{2} + 3 p_{3} t^{2} \left(t - 1\right) - p_{4} t^{3} + px\right)
Es werden jedoch nur die oben berechneten zwei Formeln verwendet, und die unklaren Punkte $ p_2 $ und $ p_3 $ werden auf 1,0 gesetzt, und jeder Term wird berechnet. Ich werde versuchen, es so C / C ++ wie möglich zu schreiben.
def lsm(points):
'''Gibt die Lösung der Methode der kleinsten Quadrate als Matrix zurück'''
p1, p4 = points[0][1], points[-1][1]
a = b = c = d = e = f = 0.
for _t, _px in points[1:-1]:
_t_1_d = (_t - 1.) ** 2
_t_1 = (_t - 1.)
_p1 = p1 * _t_1
_p2 = 3 * _t * _t_1_d
_p3 = 3 * _t ** 2 * _t_1
_p4 = p4 * _t ** 3
_j_p2 = 3 * _t * _t_1_d
_j_p3 = 3 * _t ** 2 * _t_1
a += _p2 * _j_p2
b += -_p3 * _j_p2
c += -_p2 * _j_p3
d += _p3 * _j_p3
e += -((-_p1 + _p4 - _px) * _j_p2)
f += -((_p1 - _p4 + _px) * _j_p3)
return [
[a, b],
[c, d]
], [e, f]
def inv(mat, vec):
'''Löse simultane Gleichungen mit inverser Matrix'''
dm = mat[0][0] * mat[1][1] - mat[1][0] * mat[0][1]
assert (dm != 0.0)
a = mat[1][1] / dm
b = -mat[0][1] / dm
c = -mat[1][0] / dm
d = mat[0][0] / dm
return (vec[0] * a + vec[1] * b,
vec[0] * c + vec[1] * d)
points = [(0, 0), (0.2, 65), (0.7, 45), (1.0, 100)]
inv(*lsm(points))
Nachwort
Die Methode der kleinsten Quadrate ist eine sehr leistungsfähige Technik, die angenähert werden kann, wenn es sich um eine differenzierbare Funktion handelt.
Die Annäherung der Bezier-Kurve scheint Spaß zu machen, wenn sie auf ein Spielprogramm angewendet wird. Wenn Sie die partielle Ableitung mit "Sympy" im Voraus berechnen, können Sie sie in anderen Sprachen wie C / C ++ berechnen, indem Sie einfach einen Algorithmus implementieren, der simultane Gleichungen löst. (Es gibt viele Algorithmen für simultane Gleichungen, wie zum Beispiel "Gauss Jordan Sweeping Method")
Schließlich ist "Sympy" ein sehr leistungsfähiges Werkzeug, das Formeln mit vielen Variablen leicht optimiert / berechnet. Ich kann nicht gut rechnen, deshalb ist es sehr nützlich.