Approximer une courbe de Bézier à travers un point spécifié en utilisant la méthode des moindres carrés en Python
Le déclencheur était que lorsque j'organisais les programmes que j'avais écrits dans le passé, l'ajustement de la courbe de Bézier (Curve fit % 9A% E3% 81% 82% E3% 81% A6% E3% 81% AF% E3% 82% 81)) est sorti, et je pensais que c'était nostalgique, mais la formule était difficile. Je ne pouvais pas du tout comprendre.
Pourquoi puis-je l'approcher avec une formule aussi mystérieuse, je suis incroyable dans l'ancien temps.
Apparemment, il a été approximé en utilisant la méthode des moindres carrés, je voudrais l'expliquer en Python après avoir examiné. J'ai sorti le livre de référence et l'ai réétudié. Au fait, la référence est du Dr. Kenichi Kanaya ["Des mathématiques optimisées que vous pouvez comprendre - des principes de base aux méthodes de calcul"](https://www.amazon.co.jp/ % 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). Je recommande ce livre car il est très facile à comprendre.
Quelle est la méthode des moindres carrés?
Méthode du carré minimum C'est une méthode typique pour approximer des données et des fonctions complexes, et c'est la base la plus importante pour l'analyse des données. Le champ d'application est très large et c'est une méthode à la fois pratique et esthétique que toute fonction peut être approchée si elle peut être différenciée.
Par exemple, disons que vous avez 20 points ici.
Mettons une ligne droite qui passe au plus près des 20 points.
Pour la fonction linéaire $ f (x) = ax + b $, trouvez les paramètres $ a $, $ b $ qui passent le plus proche de tous les points. La méthode typique pour trouver $ a $ et $ b $ est la ** méthode des moindres carrés **.
On l'appelle la méthode des moindres carrés car elle se rapproche de l'emplacement spécifié et de l'emplacement tracé de sorte que la somme des carrés des ** erreurs soit minimisée **.
Approximation en ligne droite
Tout d'abord, essayons l'approximation de ligne droite la plus simple et dessinons une ligne droite qui se rapproche des quatre points suivants.
points = [(107, 101), (449, 617), (816, 876), (1105,1153)]
La fonction linéaire est $ f (x) = ax + b $, $ a $ est la pente et $ b $ est la section. Pour la fonction $ f (x) $, trouvez les paramètres $ a $ et $ b $ qui minimisent le carré de l'erreur au point spécifié.
En supposant que l'emplacement spécifié est $ y $, la fonction $ j (x, y) $ qui calcule le carré de l'erreur est la suivante.
j(x, y) = (y - (ax + b))^2
Puisque la somme de $ j $ doit être minimisée pour N points, la formule de la méthode des moindres carrés est la suivante.
J = \frac {1}{2}\sum_{k=1}^{n} j(x_k, y_k)
Trouvez la dérivée partielle de la fonction d'erreur j (x, y)
Soudain, le mot ésotérique ** fonction de biais ** est sorti, mais ne vous inquiétez de rien. Même les élèves du primaire peuvent calculer s'ils connaissent la méthode de calcul.
La différenciation d'une fonction qui a deux variables ou plus par rapport à une variable est appelée ** différenciation partielle **. Et la fonction obtenue par différenciation partielle est appelée ** dérivée partielle **.
</ i> Explication pour ceux qui ont oublié la différenciation
La différenciation est $ (X ^ n) '= nX ^ {n-1} $, et l'exposant est n comme coefficient et l'exposant est n-1.
Exemple:
La différenciation est la pente de la ligne tangente du graphique, qui peut être reformulée en trouvant la vitesse instantanée du graphique. Un point avec un différentiel de 0 est un point tournant dans le cas d'un graphe convexe.
Trouvez le différentiel partiel à l'aide de Sympy
Si vous comprenez la différenciation, il n'y a rien de nouveau dans la différenciation partielle, mais ce n'est pas grave si vous ne comprenez pas les choses difficiles. Nous avons Sympy (et Jupyter)!
Si vous ne connaissez pas Sympy, veuillez lire Créer l'environnement de calcul le plus puissant avec Sympy + Jupyter.
Commençons Jupyter tout de suite.
from sympy import *
init_session()
Ceci termine la configuration de Sympy et Jupyter. Définissez une fonction $ j $ pour trouver l'erreur.
a, b = symbols("a b")
j = (y - (a * x + b)) ** 2
Vous pouvez trouver la différenciation partielle avec la fonction sympy.diff. Premièrement, la différenciation partielle est effectuée pour $ a $.
La fonction de biais pour $ a $ s'écrit $ \ frac {\ partial j} {\ partial a} $, et pour $ b $, elle s'écrit $ \ frac {\ partial j} {\ partial b} $.
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
J'ai pu trouver la dérivée partielle de la fonction linéaire, et avec Sympy, ce serait une victoire facile.
La méthode du carré minimum correspond à 1/2 de l'erreur totale, calculons-la donc.
L'attribution d'une valeur à une expression se fait avec la fonction 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
Vous avez 1/2 de chaque somme, et tout ce que vous avez à faire est de résoudre les équations simultanées de ces deux équations.
\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}
La réponse est «a = 1.01694642025091», «b = 57.0059292596265». Remplacez-le dans l'équation linéaire $ y = ax + b $. Tracons-le tout de suite.
Puisqu'il s'agit d'une ligne droite, elle ne peut pas passer par tous les points, mais elle passe par un endroit proche d'elle.
Approximation de la courbe de Bézier
Nous avons approché une ligne droite simple plus tôt, mais la courbe de Bézier d'une courbe cubique peut être approchée exactement par la même procédure.
Dans la courbe de Bézier, $ p_1 $ est le point de départ, $ p_4 $ est le point final, et $ p_2 $ et $ p_3 $ sont les points de contrôle, respectivement.
La formule de la courbe de Bézier cubique est:
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}
La plage de t est 0 ≤ t ≤ 1.
En supposant que la formule de la courbe de Bézier est $ bz $ et que le point à approximer est $ p_x $, la fonction $ j $ pour calculer l'erreur est la suivante.
j = (p_x - bz)^{2}
La définition de la méthode des moindres carrés était la suivante.
J = \frac {1}{2}\sum_{k=1}^{n} j(x_k, y_k)
Les paramètres que vous voulez trouver sont $ p_2 $ et $ p_3 $, alors faites une différenciation partielle avec chaque paramètre. Il y a 1/2 de la définition de la méthode des moindres carrés, donc multiplions la fonction de déviation par 1/2 à l'avance.
Écrivons avec Sympy à la fois.
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)
1/2 des fonctions de biais de $ p_2 $ et $ p_3 $ sont les suivantes, respectivement.
\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)
Définissez la valeur de paire du point $ p_x $ que vous voulez passer et $ t $ à ce moment. Par exemple, supposons que vous ayez les 4 points suivants.
points = [(0, 0), (0.2, 65), (0.7, 45), (1.0, 100)]
Les points de début et de fin sont $ p_1 $ et $ p_4 $, donc ces deux valeurs sont connues. Pour $ p_2 $ et $ p_3 $ que vous voulez trouver, résolvez la somme des dérivées partielles avec une équation de série.
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}
La réponse est «p2 = 180,456349206349, p3 = −53,0753968253968».
Vous pouvez l'approcher correctement!
Regardons divers paramètres.
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)]
À moins que ce ne soit une valeur très extrême, il semble que cela puisse passer correctement. Bien sûr, il y a certaines valeurs qui ne peuvent pas être transmises, mais d'une manière ou d'une autre, cela trace une ligne comme celle-là.
Essayez de générer diverses courbes
Enfin, sur Jupyter, la courbe de Bézier est approchée par la méthode du carré minimum, et les cellules à tracer sont décrites.
%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])
Si vous entrez une liste de points par lesquels la fonction plot doit passer, elle fera de son mieux pour passer à travers ces points en utilisant la méthode des moindres carrés.
L'utilisation est appelée comme suit.
points = [(0, 0), (0.1, -40), (0.5, -50), (0.8, 30), (1.0, 100)]
plot(points)
Résoudre sans utiliser Sympy
Ecrivez une implémentation dans un environnement sans Sympy tel que 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)
Cependant, seules les deux formules calculées ci-dessus sont utilisées, et les points peu clairs $ p_2 $ et $ p_3 $ sont définis comme 1.0, et chaque terme est calculé. J'essaierai de l'écrire en C / C ++ que possible.
def lsm(points):
'''Renvoie la solution de la méthode des moindres carrés sous forme de matrice'''
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):
'''Résoudre des équations simultanées avec matrice inverse'''
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))
Épilogue
La méthode des moindres carrés est une technique très puissante qui peut être approchée s'il s'agit d'une fonction différenciable.
L'approximation de la courbe de Bézier semble amusante lorsqu'elle est appliquée à un programme de jeu.
Si vous calculez à l'avance la dérivée partielle avec Sympy, vous pouvez la calculer dans d'autres langages tels que C / C ++ en implémentant simplement un algorithme qui résout les équations simultanées.
(Il existe de nombreux algorithmes pour les équations simultanées, comme "Gauss Jordan Sweeping Method")
Après tout, Sympy est un outil très puissant, et il optimise / calcule facilement des formules avec de nombreuses variables.
Je ne suis pas doué pour le calcul, donc c'est très utile.
Recommended Posts