Aproximación con polinomios

Buscamos un polinomio de grado a lo sumo $n$: $$ p(x) = a_0 + a_1\cdot x + \dots + a_n\cdot x^n $$ pero tenemos $m+1>n+1$ puntos de interpolación $$ (x_0,y_0),\:(x_1,y_1),\dots,(x_m,y_m). $$ Sabemos que en general son demasiadas condiciones, y no hay ningún polinomio que pase por todos los puntos (aunque puede ocurrir en casos concretos).

Es decir, ningún polinomio cumple estas $m+1$ condiciones: $$ p(x_0)=y_0,\:p(x_1)=y_1,\dots,p(x_m)=y_m. $$

  • Aproximación: dataset: $\{(x_i,y_i)\}_{i=0}^m$ buscamos una función que no pasa exactamente por los n puntos, y tenemos menos parámetros que datos

El objetivo al hacer aproximación (con polinomios) es minimizar el error total, usualmente el error cuadrático total:

$$ \sum_{i=0}^m \left(y_i - P(x_i)\right)^2 $$

para $$ P(x)=\sum_{j=0}^n a_j (x)^j $$

  • ¡¡Es una función cuadrática (y convexa) de las incógnitas $a_j$!! => podemos encontrar el mínimo global fácilmente con técnicas de Cálculo II, que nos da la "solución de mínimos cuadrados", que también estudiasteis en Álgebra Lineal.

Ejemplo de ajuste lineal

Tratamos de predecir la temperatura a la que hierve el agua ( T ), conocida la presión atmosférica ( P ) en el lugar y momento en que hacemos el experimento.

Para ello, contamos con un conjunto de mediciones de ambas variables, con la temperatura en grados Fahrenheit y la presión en pulgadas de mercurio (sea lo que sea, es una unidad de medidad de presión). Por ejemplo, en un cierto punto de los Alpes, un cierto día, el barómetro marcaba 20.79 pulgadas de mercurio, y el agua hirvió a 194.5 grados Fahrenheit. Las mediciones se realizaron en el mismo lugar geográfico, pero en días distintos, con distintas condiciones atmosféricas y quizá incluso por personas distintas. En estas condiciones, es imposible que ningún modelo prediga con exactitud el valor de T en función de P, pero esperamos que lo haga con un margen de error moderado.

Los datos parecen estar dispuestos aproximadamente sobre una recta, de modo que intentamos ajustar un modelo lineal, de la forma:

$$ T=a+b\cdot P $$

Nuestras incógnitas son los valores de a y b. Si el modelo lineal fuera exacto, se cumplirían las $n$ ecuaciones $$ T_i=a+b\cdot P_i,\quad i=1,\dots,n $$ donde $n=17$ es el número de mediciones. Se trata de un sistema de 17 ecuaciones para dos incógnitas, que no tiene solución. Para unos valores fijos de a y b , cometeremos un error en la medición j-ésima que será exactamente: $|T_j-(a+b\cdot P_j)|$.

Podemos escoger los valores de a y b para los que el error máximo cometido es menor, o aquellos para los que el error medio cometido es menor, o según otros muchos criterios. Es bastante habitual en estadística buscar los valores de a y b que hacen mínimo el error cuadrático total :

$$ E=\sum_j (T_j-(a+b\cdot P_j))^2 $$

In [1]:
T_values = [194.5, 194.3, 197.9, 198.4, 199.4, 199.9, 200.9, 201.1, 201.4, 201.3, 203.6, 204.6, 209.5, 208.6, 210.7, 211.9, 212.2];
P_values = [20.79, 20.79, 22.4, 22.67, 23.15, 23.35, 23.89, 23.99, 24.02, 24.01, 25.14, 26.57, 28.49, 27.76, 29.04, 29.88, 30.06];

plot(P_values, T_values,'r.');

xlabel('Presion (mm Hg)')
ylabel('Temperatura (F)')

Ejercicios

  • Investiga la documentación de polyfit para calcular el polinomio aproximador.
  • Usa polyfit para ajustar un polinomio de grado 1, y otro de grado 2, a los datos anteriores.
In [ ]:
In [ ]:
In [ ]:

Ejercicio 2

Repetimos ahora el ejercicio anterior para una serie de experimentos en un canal de ensayos en las que se mide la resistencia de varias embarcaciones similares a distinto número de Froude.

El objetivo es predecir la resistencia como función polinómica del número de Froude: $$ R = P(F) $$

Resistencia y número de Froude en las pruebas

Froude   resistence
0.125    0.03
0.15     0.18
0.175    0.4
0.2      0.73
0.225    1.3
0.25     2.16
0.275    3.35
0.3      5.06
0.325    7.14
0.35     10.36
0.375    15.25
0.4      23.15
0.425    34.62
0.45     51.5

Ejercicio

  • Usa polyfit para ajustar polinomios de grado 1, 2, 3 y 4 a los datos anteriores.
In [ ]:
In [ ]:

Construir el polinomio aproximador

Planteamos un sistema de ecuaciones para los coeficientes de un polinomio: $$ P(x) = a_0 + a_1\cdot x + \dots + a_n\cdot x^n $$ que corresponden a "el polinomio pasa por los $m+1$ puntos de interpolación" $$ P(x_0)=y_0,\:P(x_1)=y_1,\dots,P(x_m)=y_m. $$ Desarrollando esta expresión, los coeficientes del polinomio interpolador (que en general no existe si $m>n$) son la solución del sistema de ecuaciones: $$ V\cdot \mathbf{a} = \mathbf{y}, $$ donde $$ {\displaystyle V={\begin{bmatrix}1&x _{0}&x _{0}^{2}&\dots &x _{0}^{n}\\1&x _{1}&x _{1}^{2}&\dots &x _{1}^{n}\\1&x _{2}&x _{2}^{2}&\dots &x _{2}^{n}\\\vdots &\vdots &\vdots &\ddots &\vdots \\1&x _{m}&x _{m}^{2}&\dots &x _{m}^{n}\end{bmatrix}}.} $$ $$ {\displaystyle \mathbf{a}={\begin{bmatrix}a _{0}\\a_1\\\vdots\\a_n\end{bmatrix}}; \mathbf{y}={\begin{bmatrix}y_{0}\\y_1\\\vdots\\y_m\end{bmatrix}},} $$

Solución de mínimos cuadrados

El sistema de ecuaciones planteado en general no tiene solución

$$ V\cdot \mathbf{a} = \mathbf{y}, $$

Sin embargo, la solución de mínimos cuadrados está definido por

$$ \left(V^T\cdot V\right)\cdot \mathbf{a} = V^T\cdot \mathbf{y}, $$

y este sistema de ecuaciones sí tiene solución única

Ejercicio

  • Replica la solución obtenida mediante polyfit en alguno de los ejercicios anteriores planteando y resolviendo este sistema de ecuaciones.
In [ ]:
In [ ]: