Ajustement de polynômes à des données

Question

Est-il possible, étant donné un ensemble de valeurs (x,f(x)), pour trouver le polynôme de degré qui correspond le mieux aux données?

Je sais que l'interpolation polynomiale, ce qui est pour trouver un polynôme de degré n accordée n+1 de points de données, mais il y a ici un grand nombre de valeurs et nous voulons trouver un faible degré de polynôme (trouver le meilleur ajustement linéaire, le meilleur quadratique, le meilleur cubes, etc.). Elle pourrait être liée à de moindres carrés...

Plus généralement, je voudrais savoir la réponse, lorsque nous avons un multivariée de la fonction -- des points tels que l' (x,y,f(x,y)), le dire-et voulez trouver le meilleur polynôme (p(x,y)) d'un degré donné dans les variables. (Plus précisément un polynôme, pas des splines ou des séries de Fourier.)

À la fois la théorie et le code/les bibliothèques (de préférence en Python, mais toute langue est d'accord) serait utile.

Answer 1

Merci pour tous les réponses. Voici une autre tentative de les résumer. Pardon si je dis trop "évident" les choses: je ne savais rien à propos de la méthode des moindres carrés avant, tout était nouveau pour moi.

PAS polynôme d'interpolation

L'interpolation polynomiale est le montage d'un polynôme de degré n accordée n+1 de points de données, p. ex. la recherche d'un cube qui passe exactement par quatre points donnés. Comme dit dans la question, ce n'était pas ce que je voulais-j'ai eu beaucoup de points et je voulais un petit polynomiale de degré (qui ne sera d'environ ajustement, sauf si nous avons eu de la chance)-mais depuis quelques réponses a insisté sur le fait d'en parler, je dois mentionner :) polynôme de Lagrange, de la matrice de Vandermonde, etc.

Qu'est-ce que la méthode des moindres carrés?

"Des moindres carrés" est une définition particulière/critère/"métrique" de "comment" d'un polynôme s'adapte. (Il y en a d'autres, mais c'est plus simple.) Disons que vous essayez d'ajustement d'un polynôme p(x,y) = a + bx + cy + dx² + ey² + fxy à certains points de données (x_i,y_i,Z_i) (où "Z_i" a "f(x_i,y_i)" dans la question). Avec la méthode des moindres carrés le problème est de trouver le "meilleur" des coefficients (a,b,c,d,e,f), tel que ce qui est réduite gardé "moins") est la "somme des carrés des résidus", à savoir

S = ∑_i (a + bx_i + cy_i + dx_i² + ey_i² + fx_iy_i - Z_i)²

La théorie de l'

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

∑_i2(a + ... + fx_iy_i - Z_i)x_iy_i = 0

et similaires équations pour a, b, c, d, e. Notez que ce sont juste des équations linéaires à une...f. Donc, nous pouvons les résoudre avec l'élimination de Gauss ou toute les méthodes habituelles.

C'est encore appelé "linéaire des moindres carrés", parce que, bien que la fonction que nous voulions, c'était un polynôme quadratique, il est toujours linéaire dans les paramètres (a,b,c,d,e,f). Notez que la même chose travaille quand on veut p(x,y) toute "combinaison linéaire" de l' arbitraire des fonctions f_j, au lieu de simplement un polynôme (= "combinaison linéaire de monomials").

Code

Pour le univariée cas (quand il y a seulement de la variable x - f_j sont monomials x^j), il est Numpy l' polyfit:

>>> 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 multivariée cas, ou linéaire des moindres carrés en général, il est SciPy. Comme expliqué dans la documentation, il faut une matrice A de la valeurs de f_j(x_i). (La théorie est qu'il trouve le Moore-Penrose pseudo-inverse de A.) Avec notre exemple ci-dessus impliquant (x_i,y_i,Z_i), l'ajustement d'un polynôme signifie que la f_j sont les monomials x⁽⁾y⁽⁾. La suite de trouver la meilleure quadratique (ou polynôme de tout autre diplôme, si vous changez le "degré = 2" de la ligne):

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

imprime

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

ainsi, il a découvert que le polynôme x²+2xy+y²+0.01. [Le dernier terme est parfois -0.01 et parfois 0, ce qui est prévisible à cause du bruit aléatoire, nous avons ajouté.]

Des Alternatives à Python+Numpy/Scipy sont R et les Systèmes de calcul formel: la Sauge, Mathematica, Matlab, Maple. Excel peut être en mesure de le faire. Numérique Recettes décrit les méthodes à mettre en œuvre nous-mêmes (en C, Fortran).

Préoccupations

• Elle est fortement influencée par la façon dont les points sont choisis. Quand j'ai eu x=y=range(20) au lieu des points aléatoires, il a toujours produit de 1,33 x²+1.33 xy+1.33 y², ce qui est déroutant... jusqu'à ce que j'ai réalisé que parce que j'ai toujours eu x[i]=y[i], les polynômes sont les mêmes: x²+2xy+y² = 4x² = (4/3)(x²+xy+y²). Donc la morale est qu'il est important de choisir les points avec soin pour obtenir le "droit" polynôme. (Si vous pouvez choisir, vous devez choisir de Tchebychev nœuds pour l'interpolation polynomiale; vous ne savez pas si la même chose est vraie pour le moins de carrés ainsi.)
  
• Le surajustement: la hausse des polynômes de degré peut toujours ajuster les données mieux. Si vous modifiez l' degree à 3, 4 ou 5, il est encore principalement reconnaît le même polynôme quadratique (les coefficients sont 0 pour de plus haut degré termes), mais pour les plus grands degrés, il commence à montage supérieur, les polynômes de degré. Mais même avec le degré 6, de prendre de plus grandes n (plus de points de données au lieu de 20, de dire 200) s'inscrit toujours le polynôme quadratique. Donc la morale est pour éviter le surajustement, pour lesquels il peut être utile de prendre autant de points de données que possible.
• Il pourrait y avoir des problèmes de stabilité numérique je ne comprends pas tout.
• Si vous n'avez pas besoin d'un polynôme, vous pouvez obtenir un meilleur cadre avec d'autres types de fonctions, par exemple, les splines (polynômes par morceaux).

Answer 2

Oui, la façon dont cela est fait habituellement à l'est par l'utilisation de la méthode des moindres carrés. Il y a d'autres façons de spécifier comment un polynôme convient, mais que la théorie de la solution la plus simple est la méthode des moindres carrés. La théorie générale est appelée régression linéaire.

Votre meilleur pari est probablement commencer avec Recettes ou Numérique.

R est gratuit, et de faire ce que vous voulez et plus, mais il a une courbe d'apprentissage importante.

Si vous avez accès à Mathematica, vous pouvez utiliser la fonction Fit faire un ajustement des moindres carrés. J'imagine que Matlab et de son open source homologue Octave ont une fonction similaire.

Answer 3

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)
 

Answer 4

Si vous souhaitez adapter l' (xi, f(xi)) pour un polynôme de degré n , alors vous mettre en place un linéaire des moindres carrés problème avec les données (1, xi, xi, xi^2, ..., xi^n, f(xi) ). Cela renvoie un ensemble de coefficients (c0, c1, ..., cn) de sorte que la meilleure polynôme est *y = c0 + c1 * x + c2 * x^2 + ... + cn * x^n.*

Vous pouvez généraliser ce soit deux de plus que d'une variable dépendante y compris par les puissances de y et les combinaisons de x et de y dans le problème.

Answer 5

Les polynômes de Lagrange (comme @j w affichés) de vous donner un ajustement exact sur les points que vous indiquez, mais avec des polynômes de degré plus que de 5 ou 6, vous pouvez exécuter en instabilité numérique.

Moins de carrés vous donne le meilleur ajustement polynomial avec l'erreur définie comme la somme des carrés des erreurs individuelles. (prendre de la distance le long de l'axe entre les points que vous avez et la fonction que les résultats, les carrés entre eux, et de les résumer) Le MATLAB polyfit fonction de cela, et avec de multiples retourner les arguments, vous pouvez automatiquement prendre soin de mise à l'échelle/compenser les problèmes (par exemple, si vous avez 100 points entre x=312.1 et 312.3, et que vous voulez une 6e polynomiale de degré, vous allez avoir à calculer u = (x-312.2)/0.1 donc les valeurs u sont distribués entre -1 et +=).

NOTEZ que les résultats de la méthode des moindres carrés ajustements sont fortement influencés par la distribution des valeurs d'axe des abscisses. Si les valeurs de x sont également espacés, puis vous obtenez des erreurs importantes sur les extrémités. Si vous avez un cas où vous pouvez choisir les valeurs de x et vous vous souciez de la déviation maximale de votre fonction connue et une interpolation polynomiale, alors l'utilisation de polynômes de Tchebychev va vous donner quelque chose qui est proche du parfait polynôme minimax (ce qui est très difficile à calculer). Cette question est examinée en détail dans les Recettes ou Numérique.

Edit: De ce que je comprends, tout cela fonctionne bien pour les fonctions d'une seule variable. Pour les fonctions à plusieurs variables, il est susceptible d'être beaucoup plus difficile si le diplôme est de plus de, disons, 2. J'ai trouvé une référence sur Google Livres.