Adapter des polynômes à des données

Y a-t-il un moyen, étant donné un ensemble de valeurs (x,f(x)) , de trouver le polynôme d’un degré donné qui correspond le mieux aux données?

Je connais l’ interpolation polynomiale , qui sert à trouver un polynôme de degré n donné n+1 points de données, mais ici il y a un grand nombre de valeurs et nous voulons trouver un polynôme de faible degré (trouver le meilleur ajustement linéaire, le meilleur cubique, etc.). Cela peut être lié aux moindres carrés …

Plus généralement, je voudrais connaître la réponse lorsque nous avons une fonction multivariée – des points comme (x,y,f(x,y)) , par exemple – et que nous voulons trouver le meilleur polynôme ( p(x,y) ) d’un degré donné dans les variables. (Spécifiquement un polynôme, pas des splines ou des séries de Fourier.)

La théorie et le code / bibliothèques (de préférence en Python, mais n’importe quel langage est correct) seraient utiles.

Merci pour les réponses de tous. Voici une autre tentative de les résumer. Pardon si je dis trop de choses “évidentes”: je ne savais rien des moindres carrés avant, alors tout était nouveau pour moi.

PAS l’interpolation polynomiale

L’interpolation polynomiale adapte un polynôme de degré n n+1 points de données, par exemple en trouvant un cube qui passe exactement à travers quatre points donnés. Comme je le disais dans la question, ce n’était pas ce que je voulais – j’avais beaucoup de points et je voulais un polynôme de petit degré (qui ne conviendrait qu’à peu près , sauf si nous avons eu de la chance). à ce sujet, je devrais les mentionner 🙂 Polynôme de Lagrange , masortingce de Vandermonde , etc.

Qu’est-ce que les moindres carrés?

Les “moindres carrés” sont une définition / un critère / une “mésortingque” particulière de “comment bien” correspond un polynôme. (Il y en a d’autres, mais c’est le plus simple.) Disons que vous essayez d’adapter un polynôme p (x, y) = a + bx + cy + dx 2 + ey 2 + fxy à des points de données donnés (x i , y i , Z i ) (où “Z i ” était “f (x i , y i )” dans la question). Avec les moindres carrés, le problème est de trouver les “meilleurs” coefficients (a, b, c, d, e, f), de sorte que ce qui est minimisé (conservé “le moins”) soit la “sum des résidus au carré”, à savoir

S = ∑ i (a + bx i + cy i + dx i 2 + ey i 2 + fx i y i – Z i ) 2

Théorie

L’idée importante est que si vous regardez S en fonction de (a, b, c, d, e, f), alors S est minimisé à un point où son gradient est 0 . Cela signifie que, par exemple, ∂S / ∂f = 0, c’est-à-dire que

2 i 2 (a +… + fx i y i – Z i ) x i y i = 0

et équations similaires pour a, b, c, d, e. Notez que ce ne sont que des équations linéaires dans un… f. Nous pouvons donc les résoudre avec l’élimination gaussienne ou l’une des méthodes habituelles .

Ceci est encore appelé “moindres carrés linéaires”, car bien que la fonction que nous voulions soit un polynôme quadratique, elle est toujours linéaire dans les parameters (a, b, c, d, e, f). Notez que la même chose fonctionne lorsque l’on veut que p (x, y) soit une “combinaison linéaire” de fonctions arbitraires f j , au lieu d’un simple polynôme (= “combinaison linéaire de monômes”).

Code

Pour le cas univarié (quand il n’y a que la variable x – les f j sont des monômes x j ), il y a polyfit de Numpy:

 >>> import numpy >>> xs = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] >>> ys = [1.1, 3.9, 11.2, 21.5, 34.8, 51, 70.2, 92.3, 117.4, 145.5] >>> p = numpy.poly1d(numpy.polyfit(xs, ys, deg=2)) >>> print p 2 1.517 x + 2.483 x + 0.4927 

Pour le cas multivarié, ou les moindres carrés linéaires en général, il y a SciPy. Comme expliqué dans sa documentation , il prend une masortingce A des valeurs f j ( x i ). (La théorie est qu’elle trouve la pseudoinverse de Moore-Penrose de A.) Avec notre exemple ci-dessus impliquant (x i , y i , Z i ), ajuster un polynôme signifie que f j sont les monômes x () y () . Ce qui suit trouve le meilleur quadratique (ou meilleur polynôme de tout autre degré, si vous modifiez la ligne “degree = 2”):

 from scipy import linalg import random n = 20 x = [100*random.random() for i in range(n)] y = [100*random.random() for i in range(n)] Z = [(x[i]+y[i])**2 + 0.01*random.random() for i in range(n)] degree = 2 A = [] for i in range(n): A.append([]) for xd in range(degree+1): for yd in range(degree+1-xd): A[i].append((x[i]**xd)*(y[i]**yd)) #f_j(x_i) c,_,_,_ = linalg.lstsq(A,Z) j = 0 for xd in range(0,degree+1): for yd in range(0,degree+1-xd): print " + (%.2f)x^%dy^%d" % (c[j], xd, yd), j += 1 

estampes

  + (0.01)x^0y^0 + (-0.00)x^0y^1 + (1.00)x^0y^2 + (-0.00)x^1y^0 + (2.00)x^1y^1 + (1.00)x^2y^0 

il a donc découvert que le polynôme est x 2 + 2xy + y 2 + 0,01. [Le dernier terme est parfois -0.01 et parfois 0, ce à quoi on peut s’attendre en raison du bruit aléatoire que nous avons ajouté.]

Les alternatives à Python + Numpy / Scipy sont les systèmes algébriques R et Computer: Sage , Mathematica, Matlab, Maple. Même Excel peut le faire. Recettes numériques discute des méthodes pour l’implémenter nous-mêmes (en C, Fortran).

Les soucis

  • Il est fortement influencé par la manière dont les points sont choisis . Quand j’avais x=y=range(20) au lieu des points aléatoires, il produisait toujours 1,33x 2 + 1.33xy + 1.33y 2 , ce qui était déroutant … jusqu’à ce que je m’en rende compte parce que j’avais toujours x[i]=y[i] , les polynômes étaient les mêmes: x 2 + 2xy + y 2 = 4x 2 = (4/3) (x 2 + xy + y 2 ). La morale est donc qu’il est important de choisir les points avec soin pour obtenir le “bon” polynôme. (Si vous pouvez choisir, vous devriez choisir les nœuds Chebyshev pour l’interpolation polynomiale; vous ne savez pas si la même chose est vraie pour les moindres carrés.)
  • Surajustement : les polynômes de haut degré peuvent toujours mieux adapter les données. Si vous modifiez le degree à 3 ou 4 ou 5, il reconnaît toujours le même polynôme quadratique (les coefficients sont égaux à 0 pour les termes plus élevés), mais pour les degrés supérieurs, il commence à ajuster les polynômes de degré supérieur. Mais même avec le degré 6, prendre plus grand n (plus de points de données au lieu de 20, disons 200) correspond toujours au polynôme quadratique. La morale consiste donc à éviter le sur-ajustement, pour lequel il pourrait être utile de prendre autant de points de données que possible.
  • Il pourrait y avoir des problèmes de stabilité numérique que je ne comprends pas bien.
  • Si vous n’avez pas besoin d’un polynôme, vous pouvez obtenir de meilleurs ajustements avec d’autres types de fonctions, par exemple les splines (polynômes par morceaux).

Oui, la manière la plus courante consiste à utiliser les moindres carrés. Il y a d’autres manières de spécifier le bon ajustement d’un polynôme, mais la théorie est la plus simple pour les moindres carrés. La théorie générale est appelée régression linéaire.

Votre meilleur pari est probablement de commencer avec des recettes numériques .

R est gratuit et fera tout ce que vous voulez et plus, mais il a une grande courbe d’apprentissage.

Si vous avez access à Mathematica, vous pouvez utiliser la fonction Fit pour effectuer un ajustement par moindres carrés. J’imagine que Matlab et son homologue open source Octave ont une fonction similaire.

Pour (x, f (x)) cas:

 import numpy x = numpy.arange(10) y = x**2 coeffs = numpy.polyfit(x, y, deg=2) poly = numpy.poly1d(coeffs) print poly yp = numpy.polyval(poly, x) print (yp-y) 

N’oubliez pas qu’un polynôme de degré supérieur correspond TOUJOURS aux données. Les polynômes de degré plus élevé conduisent généralement à des fonctions hautement improbables (voir le razor d’Occam ), bien que sur-ajusté. Vous voulez trouver un équilibre entre simplicité (degré de polynôme) et ajustement (par exemple, erreur de moindre carré). Quantitativement, il existe des tests pour cela, le critère d’information d’Akaike ou le critère d’information bayésien . Ces tests donnent une note sur le modèle à privilégier.

Si vous voulez adapter le (xi, f (xi)) à un polynôme de degré n, vous devez définir un problème de moindres carrés linéaire avec les données (1, xi, xi, xi ^ 2, …, xi ^ n, f (xi)). Cela renverra un ensemble de coefficients (c0, c1, …, cn) pour que le meilleur polynôme soit * y = c0 + c1 * x + c2 * x ^ 2 + … + cn * x ^ n. *

Vous pouvez généraliser ces deux variables dépendantes en incluant des puissances de y et des combinaisons de x et y dans le problème.

Les polynômes de Lagrange (tels que @jw posted) vous permettent d’ajuster exactement les points que vous spécifiez, mais avec des polynômes de degré supérieur à 5 ou 6, vous pouvez rencontrer une instabilité numérique.

Les moindres carrés vous donnent le polynôme «meilleur ajustement» avec une erreur définie comme la sum des carrés des erreurs individuelles. (prenez la distance le long de l’axe y entre les points que vous avez et la fonction qui en résulte, polyfit les et additionnez-les) La fonction polyfit MATLAB le fait, et avec plusieurs arguments de retour, / offset (par exemple, si vous avez 100 points entre x = 312.1 et 312.3, et que vous voulez un polynôme du 6ème degré, vous voudrez calculer u = (x-312.2) /0.1 pour que les valeurs-u soient dissortingbuées entre -1 et + =).

Notez que les résultats des ajustements par les moindres carrés sont fortement influencés par la dissortingbution des valeurs de l’axe des x. Si les valeurs x sont espacées de manière égale, vous obtiendrez des erreurs plus importantes aux extrémités. Si vous avez un cas où vous pouvez choisir les valeurs x et vous soucier de l’écart maximal par rapport à votre fonction connue et d’un polynôme d’interpolation, alors l’utilisation des polynômes de Chebyshev vous donnera quelque chose qui se rapproche difficile à calculer). Ceci est discuté en détail dans Recettes numériques.

Edit: D’après ce que je comprends, tout fonctionne bien pour les fonctions d’une variable. Pour les fonctions à plusieurs variables, cela risque d’être beaucoup plus difficile si le diplôme est supérieur à, disons, 2. J’ai trouvé une référence sur Google Livres .

à l’université, nous avions ce livre que je trouve toujours extrêmement utile: Conté, de Boor; parsing numérique élémentaire; Mc Grow Hill. Le paragraphe pertinent est 6.2: Montage des données.
L’exemple de code est fourni avec FORTRAN, et les listes ne sont pas très lisibles non plus, mais les explications sont profondes et claires en même temps. vous finissez par comprendre ce que vous faites, pas seulement le faire (comme c’est mon expérience des recettes numériques).
Je commence généralement par les recettes numériques, mais pour des choses comme celles-ci, je dois rapidement attraper Conte-de Boor.

peut-être mieux poster un code … c’est un peu dépouillé, mais les parties les plus pertinentes sont là. ça repose sur numpy, évidemment!

 def Tn(n, x): if n==0: return 1.0 elif n==1: return float(x) else: return (2.0 * x * Tn(n - 1, x)) - Tn(n - 2, x) class ChebyshevFit: def __init__(self): self.Tn = Memoize(Tn) def fit(self, data, degree=None): """fit the data by a 'minimal squares' linear combination of chebyshev polinomials. cfr: Conte, de Boor; elementary numerical analysis; Mc Grow Hill (6.2: Data Fitting) """ if degree is None: degree = 5 data = sorted(data) self.range = start, end = (min(data)[0], max(data)[0]) self.halfwidth = (end - start) / 2.0 vec_x = [(x - start - self.halfwidth)/self.halfwidth for (x, y) in data] vec_f = [y for (x, y) in data] mat_phi = [numpy.array([self.Tn(i, x) for x in vec_x]) for i in range(degree+1)] mat_A = numpy.inner(mat_phi, mat_phi) vec_b = numpy.inner(vec_f, mat_phi) self.coefficients = numpy.linalg.solve(mat_A, vec_b) self.degree = degree def evaluate(self, x): """use Clenshaw algorithm http://en.wikipedia.org/wiki/Clenshaw_algorithm """ x = (x-self.range[0]-self.halfwidth) / self.halfwidth b_2 = float(self.coefficients[self.degree]) b_1 = 2 * x * b_2 + float(self.coefficients[self.degree - 1]) for i in range(2, self.degree): b_1, b_2 = 2.0 * x * b_1 + self.coefficients[self.degree - i] - b_2, b_1 else: b_0 = x*b_1 + self.coefficients[0] - b_2 return b_0 

Rappelez-vous qu’il y a une grande différence entre l’ approximation du polynôme et la recherche exacte .

Par exemple, si je vous donne 4 points, vous pourriez

  1. Approximer une ligne avec une méthode comme les moindres carrés
  2. Approximer une parabole avec une méthode comme les moindres carrés
  3. Trouvez une fonction cubique exacte à travers ces quatre points.

Assurez-vous de sélectionner la méthode qui vous convient!

Il est plutôt facile de faire peur rapidement en utilisant les fonctions de masortingce d’Excel si vous savez comment représenter le problème des moindres carrés en tant que problème d’algèbre linéaire. (Cela dépend de la fiabilité avec laquelle Excel est considéré comme un résolveur d’algèbre linéaire.)

Le polynôme de lagrange est en quelque sorte le polynôme d’interpolation “le plus simple” qui correspond à un ensemble donné de points de données.

Cela pose parfois problème car il peut varier énormément entre les points de données.