math - polyfit - ajuste polinomial por minimos cuadrados matlab

Ajuste de polinomios a los datos

Question

¿Hay alguna manera, dada una serie de valores (x,f(x)) , de encontrar el polinomio de un grado dado que mejor se ajuste a los datos?

Conozco la interpolación polinómica , que es para encontrar un polinomio de grado n dado n+1 puntos de datos, pero aquí hay una gran cantidad de valores y queremos encontrar un polinomio de bajo grado (encontrar el mejor ajuste lineal, el mejor cuadrático, el mejor cúbica, etc.). Podría estar relacionado con mínimos cuadrados ...

En términos más generales, me gustaría saber la respuesta cuando tenemos una función multivariante: puntos como (x,y,f(x,y)) , por ejemplo, y queremos encontrar el mejor polinomio ( p(x,y) ) de un grado dado en las variables. (Específicamente un polinomio, no splines o series de Fourier.)

Tanto la teoría como el código / las bibliotecas (preferiblemente en Python, pero cualquier lenguaje está bien) serían útiles.

Answer 1

El polinomio de Lagrange es en cierto sentido el polinomio de interpolación "más simple" que se ajusta a un conjunto dado de puntos de datos.

A veces es problemático porque puede variar enormemente entre los puntos de datos.

Answer 2

Es bastante fácil asustar a un ajuste rápido utilizando las funciones de matriz de Excel si sabe cómo representar el problema de mínimos cuadrados como un problema de álgebra lineal. (Eso depende de cuán confiable piense que Excel es como un solucionador de álgebra lineal.)

Answer 3

Gracias por las respuestas de todos. Aquí hay otro intento de resumirlos. Perdonen si digo demasiadas cosas "obvias": antes no sabía nada sobre los mínimos cuadrados, así que todo era nuevo para mí.

NO interpolación polinómica

La interpolación polinómica se ajusta a un polinomio de grado n dado n+1 puntos de datos, por ejemplo, la búsqueda de un cubo que pasa exactamente a través de cuatro puntos determinados. Como dije en la pregunta, esto no era lo que quería, tenía muchos puntos y quería un polinomio de grado pequeño (que solo se ajustaría aproximadamente , a menos que tuviéramos suerte), pero como algunas de las respuestas insistían en hablar al respecto, debería mencionarlos :) Polinomio de Lagrange , matriz de Vandermonde , etc.

¿Qué es mínimos cuadrados?

"Mínimos cuadrados" es una definición / criterio / "medida" particular de "qué tan bien" se ajusta a un polinomio. (Hay otros, pero esto es más simple.) Supongamos que está tratando de ajustar un polinomio p (x, y) = a + bx + cy + dx ² + ey ² + fxy a algunos puntos de datos dados (x _i , y _i , Z _i ) (donde "Z _i " fue "f (x _i , y _i )" en la pregunta). Con los mínimos cuadrados el problema es encontrar los "mejores" coeficientes (a, b, c, d, e, f), de modo que lo que se minimiza (se mantiene "al menos") es la "suma de residuales cuadrados", es decir

S = Σ _i (a + bx _i + cy _i + dx _i ² + ey _i ² + fx _i y _i - Z _i ) ²

Teoría

La idea importante es que si se mira a S como una función de (a, b, c, d, e, f), entonces S se minimized en un punto en el que su gradient es 0 . Esto significa que, por ejemplo, ∂S / ∂f = 0, es decir, que

Σ _i 2 (a + ... + fx _i y _i - Z _i ) x _i y _i = 0

y ecuaciones similares para a, b, c, d, e. Tenga en cuenta que estas son solo ecuaciones lineales en a ... f. Entonces podemos resolverlos con eliminación Gaussiana o cualquiera de los métodos usuales .

Esto todavía se llama "mínimos cuadrados lineales", porque aunque la función que queríamos era un polinomio cuadrático, sigue siendo lineal en los parámetros (a, b, c, d, e, f). Tenga en cuenta que lo mismo funciona cuando queremos que p (x, y) sea cualquier "combinación lineal" de funciones arbitrarias f _j , en lugar de simplemente un polinomio (= "combinación lineal de monomios").

Código

Para el caso univariante (cuando solo hay una variable x - las f _j son monomios x ^j ), está el polyfit de 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

Para el caso multivariado, o mínimos cuadrados lineales en general, está SciPy. Como se explica en su documentación , toma una matriz A de los valores f _j ( x _i ). (La teoría es que encuentra el pseudoinverso de Moore-Penrose de A.) Con nuestro ejemplo anterior relacionado con (x _i , y _i , Z _i ), ajustar un polinomio significa que f _j son los monomios x ⁽⁾ y ⁽⁾ . Lo siguiente encuentra el mejor cuadrático (o el mejor polinomio de cualquier otro grado, si cambia la línea "grado = 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

huellas dactilares

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

por lo que ha descubierto que el polinomio es x ² + 2xy + y ² +0.01. [El último término es a veces -0.01 y a veces 0, lo cual es de esperar debido al ruido aleatorio que agregamos].

Las alternativas a Python + Numpy / Scipy son R y Computer Algebra Systems: Sage , Mathematica, Matlab, Maple. Incluso Excel podría ser capaz de hacerlo. Las Recetas Numéricas discuten los métodos para implementarlo nosotros mismos (en C, Fortran).

Preocupaciones

• Está fuertemente influenciado por cómo se eligen los puntos . Cuando tenía x=y=range(20) lugar de los puntos aleatorios, siempre producía 1.33x ² + 1.33xy + 1.33y ² , lo cual era desconcertante ... hasta que me di cuenta de eso porque siempre tenía x[i]=y[i] , los polinomios fueron los mismos: x ² + 2xy + y ² = 4x ² = (4/3) (x ² + xy + y ² ). Entonces, la moraleja es que es importante elegir cuidadosamente los puntos para obtener el polinomio "correcto". (Si puede elegir, debe elegir los nodos de Chebyshev para la interpolación polinómica; no está seguro de que lo mismo sea cierto para los mínimos cuadrados también).
  
• Sobreajuste : polinomios de alto grado siempre pueden ajustarse mejor a los datos. Si cambia el degree a 3 o 4 o 5, aún reconoce mayormente el mismo polinomio cuadrático (los coeficientes son 0 para términos de mayor grado) pero para grados más grandes, comienza a ajustar polinomios de mayor grado. Pero incluso con el grado 6, tomar n más grande (más puntos de datos en lugar de 20, digamos 200) todavía se ajusta al polinomio cuadrático. Por lo tanto, la moraleja es evitar el sobreajuste, por lo que podría ayudar tomar tantos puntos de datos como sea posible.
• Puede haber problemas de estabilidad numérica que no entiendo completamente.
• Si no necesita un polinomio, puede obtener mejores ajustes con otros tipos de funciones, por ejemplo, splines (polinomios por partes).

Answer 4

Los polinomios de Lagrange (como se publicó en @jw) le dan un ajuste exacto en los puntos que especifique, pero con polinomios de grado superior a, por ejemplo, 5 o 6, puede encontrarse con inestabilidad numérica.

Los cuadrados mínimos le dan el polinomio de "mejor ajuste" con error definido como la suma de cuadrados de los errores individuales. (toma la distancia a lo largo del eje y entre los puntos que tienes y la función que resulta, los cuadra y los suma) La función polyfit MATLAB hace esto, y con múltiples argumentos de retorno, puedes hacer que automáticamente se encargue de escalar / problemas de compensación (por ejemplo, si tiene 100 puntos entre x = 312.1 y 312.3, y quiere un polinomio de 6º grado, querrá calcular u = (x-312.2) /0.1 para que se distribuyan los valores u entre -1 y + =).

TEN EN CUENTA que los resultados de los ajustes por mínimos cuadrados están fuertemente influenciados por la distribución de los valores del eje x. Si los valores de x están equiespaciados, obtendrás errores más grandes en los extremos. Si tiene un caso en el que puede elegir los valores xy le preocupa la desviación máxima de su función conocida y un polinomio de interpolación, el uso de polinomios de Chebyshev le dará algo que se acerca al polinomio minimax perfecto (que es muy difícil de calcular). Esto se discute con cierto detalle en Recetas Numéricas.

Editar: por lo que veo, todo esto funciona bien para las funciones de una variable. Para funciones multivariables, es probable que sea mucho más difícil si el título es más que, digamos, 2. Encontré una referencia en Google Books .

Answer 5

Para (x, f (x)) caso:

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 6

Recuerde, hay una gran diferencia entre aproximar el polinomio y encontrar uno exacto .

Por ejemplo, si le doy 4 puntos, podría

1. Aproximar una línea con un método como mínimos cuadrados
2. Aproxime una parábola con un método como mínimos cuadrados
3. Encuentra una función cúbica exacta a través de estos cuatro puntos.

¡Asegúrese de seleccionar el método adecuado para usted!

Answer 7

Sí, la forma en que esto se hace típicamente es mediante el uso de mínimos cuadrados. Hay otras formas de especificar qué tan bien encaja un polinomio, pero la teoría es más simple para los mínimos cuadrados. La teoría general se llama regresión lineal.

Tu mejor opción es comenzar con Recetas Numéricas .

R es gratis y hará todo lo que quieras y más, pero tiene una gran curva de aprendizaje.

Si tiene acceso a Mathematica, puede usar la función Fit para hacer un ajuste de mínimos cuadrados. Me imagino que Matlab y su contraparte de código abierto Octave tienen una función similar.

Answer 8

Si desea ajustar el (xi, f (xi)) a un polinomio de grado n, entonces establecería un problema lineal de mínimos cuadrados con los datos (1, xi, xi, xi ^ 2, ..., xi ^ n, f (xi)). Esto devolverá un conjunto de coeficientes (c0, c1, ..., cn) de modo que el polinomio que mejor se ajuste sea * y = c0 + c1 * x + c2 * x ^ 2 + ... + cn * x ^ n. *

Puede generalizar estas dos más de una variable dependiente incluyendo poderes de y y combinaciones de x e y en el problema.

Answer 9

Ten en cuenta que un polinomio de mayor grado SIEMPRE se ajusta mejor a los datos. Sin embargo, los polinomios de mayor grado típicamente conducen a funciones altamente improbables (ver la Navaja de Occam ), (sobreajuste). Desea encontrar un equilibrio entre simplicidad (grado de polinomio) y ajuste (p. Ej., Error de mínimos cuadrados). Cuantitativamente, hay pruebas para esto, el Criterio de información de Akaike o el Criterio de información bayesiano . Estas pruebas dan una puntuación que modelo es preferible.

Answer 10

en la universidad teníamos este libro que todavía encuentro extremadamente útil: Conte, de Boor; análisis numérico elemental; Mc Grow Hill. El párrafo relevante es 6.2: Ajuste de datos.
el código de ejemplo viene en FORTRAN, y los listados tampoco son muy legibles, pero las explicaciones son profundas y claras al mismo tiempo. terminas entendiendo lo que estás haciendo, no solo haciéndolo (como es mi experiencia de Recetas Numéricas).
Normalmente empiezo con Recetas Numéricas, pero para cosas como esta, rápidamente tengo que agarrar Conte-de Boor.

tal vez sea mejor publicar algún código ... está un poco despojado, pero las partes más relevantes están ahí. se basa en numpy, obviamente!

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