Ejemplo sencillo de polinomios de tendencia 1D

Métodos potenciales de prospección, FCAG, 2024.

Dato

$\newcommand{\B}[1]{\mathbf{#1}}$ Generamos un dato sintético simple. Una señal sinusoidal con $\lambda=10$ km y de una amplitud unitaria arbitraria. Tomamos una serie de muestras $d_1,d_2,\cdots$ no equiespaciadas para digitalizar la señal. Estas muestras las consideramos como las observaciones. Además, para visualizar la señal original, consideramos un intervalo $\Delta x$ constante y adecuado según el criterio de Nyquist.

Podemos agrupar las observaciones en el vector columna $\mathbf{d}=\begin{bmatrix} d_1 \\ \vdots \\ d_N \end{bmatrix} \: .$ Por último, simulamos error en las observaciones sumando ruido Gaussiano aleatorio en cada entrada del vector $\B{d}$.

Generamos los datos necesarios:

# señal
n0  = 100
T   =  10 # km
x0  = np.linspace(0,1*T, n0)     # muestras equiespaciadas
y0  = np.sin(2 * np.pi * x0 / T) 

# observaciones
N    = 10                       # M = 9 ajusta todas las muestras.
x    = np.random.choice(x0, N)  # muestras aleatorias.
x    = np.unique(x)
d    = np.sin(2 * np.pi * x / T)

# ruido
noise = np.random.normal(size = N)

# observaciones ruidosas
d    = d + 0.4 * noise

Ajuste con polinomio de grado $M$

Ajustamos un polinomio de grado $M$ a las observaciones $\B{d}$. Por ejemplo, para un modelo lineal ($M=1$), cada muestra $d_j$ observada en $x_j$ es modelada por

$$a\,x_j+b = d_j.$$

Expresando esta relación en notación vectorial, tenemos

$$a\,\B{x}+b\,\B{1} = \B{d},$$

donde $\B{x}$ es el vector columna $N\times 1$ con las coordenadas de los puntos de observación en $\B{x}$ y $\B{1}$ es un vector columna $N\times 1$ con el valor $1$ en cada una de sus componentes.

Otra manera de expresar el problema es reuniendo los parámetros del modelo en $\B{m}=\begin{bmatrix} a \\ b\end{bmatrix}$ y escribiendo la expresión anterior como

$$\B{A}\B{m} = \B{d},$$

donde $\B{A}=[\B{x}\: \B{1}]=\begin{bmatrix} x_1 & 1 \\ \vdots & \vdots \\ x_N & 1\end{bmatrix}$ es una matriz $N\times 2$.

Podemos utilizar numpy.polyfit() para realizar el ajuste de un polinomio de grado deg = M y obtener el vector de coeficientes $\mathbf{m}$ en m:

x = ...  # vector con las coordenadas de las observaciones
d = ...  # vector con las observaciones
M = 1    # grado del polinomio a ajustar
m = numpy.polyfit(x, d, deg = M)  # vector de parámetros ajustados
...

Una vez obtenido el ajuste, podemos evaluar el polinomio en una cuadrícula más densa que la original. Con ello podemos interpolar el dato observado y visualizar el ajuste:

Visualización de los resultados

Visualizamos los ajustes para $M=1,2,3,9$.

Analizamos el error de ajuste según el grado del polinomio. El error se calcula como la norma Euclídea de la diferencia entre los vectores de las observaciones y los valores arrojados por el ajuste polinomial:

$$e = ||\B{d}-\B{A}\B{m}||_2 .$$

Como muestra el gráfico, para $M=9$ el error es cero. ¿Implica ello un buen ajuste de tendencia?

Predicción

Por último, veamos el carácter predictivo del modelo obtenido. Para ello, generamos nuevas observaciones $\hat{\B{d}}$ y analizamos el error en la predicción con el valor arrojado por el modelo que fue ajustado para las observaciones originales $\B{d}$. El error en la predicción es $e_P = ||\hat{\B{d}}-\hat{\B{A}}\B{m}||_2$, donde $\hat{\B{A}}$ es la matriz $\B{A}$ evaluada en las coordenadas de las nuevas observaciones.

¿Un buen ajuste de los datos observados implica que obtendremos luego una buena predicción?

Eso es todo por hoy.