Clase X

Oscilador armónico amortiguado

• Repaso de cuadrados mínimos para una recta
• Comentario de Causalidad y Machine Learning
• Generalización a cualquier función
• Estimación de errores a partir de residuos
• Oscilador armónico
• Oscilador armónico amortiguado

import ipywidgets
import matplotlib.pyplot as plt
import numpy as np
from scipy import integrate, optimize

plt.rc("font", size=16)

Repaso: cuadrados minimos para una recta

Por ahora, aprendimos cuadrados mínimos para una recta $y = a x + b$.

Medimos una serie de datos $(x_i, y_i)$... O, al menos, los simulamos:

x = np.arange(0, 30)  # Genero x equiespaciados
y = (9 / 5) * x + 32  # Genero y como una recta
y = y + np.random.normal(0, 1, size=x.size)  # Agrego ruido aleatorio gaussiano

y los gráficamos:

plt.errorbar(x, y, yerr=1, fmt=".")
plt.axhline(0, color="k")
plt.axvline(0, color="k")
plt.grid()

¡Parece que hay una relación entre $y$ y $x$!

Matemáticamente: $ y = f(x) $

En particular, parece una recta.

Matemáticamente: $ y = a x + b $

¿Cuáles será la recta real?

Imposible saberlo.

¿Cuál será la mejor recta a partir de estos datos?

Ahora sí. Podemos probar variar los parámetros a mano:

@ipywidgets.interact(a=(1.0, 3.0, 0.05), b=(20, 40, 2))
def _(a=1, b=20):
    plt.errorbar(x, y, yerr=1, fmt=".")
    plt.plot(x, a * x + b, color="r", linewidth=3)

    plt.grid()
    plt.axhline(0, color="k")
    plt.axvline(0, color="k")
    plt.xlim(-2, 32)
    plt.ylim(-5, 100)

Es díficil ver con más precisión la alineación de la recta con los datos.

¿Qué gráfico podemos hacer para evaluarlo mejor?

¡El de los residuos!

@ipywidgets.interact(a=(1.5, 2.5, 0.02), b=(25, 35, 1))
def _(a=1, b=20):
    fig, axes = plt.subplots(2, 1, figsize=(6, 8))

    for ax in axes:
        ax.axhline(0, color="k")
        ax.axvline(0, color="k")
        ax.grid()
        ax.set(xlim=(-2, 32))

    axes[0].set(ylim=(-5, 100))
    axes[1].set(ylim=(-5, 5))

    y_recta = a * x + b

    axes[0].errorbar(x, y, yerr=1, fmt=".")
    axes[0].plot(x, y_recta, color="r", linewidth=3)

    axes[1].errorbar(x, y - y_recta, yerr=1, fmt=".")
    axes[1].plot(x, y_recta - y_recta, "--r")

Para encontrar una solución matemática, tenemos que definir que es "mejor".

Hay distintas maneras, pero la que usamos se basa en medir la distancia cuadrática de los $y_i$ a la recta en $x_i$:

$$ \Big (y_i - (ax_i + b) \Big)^2 = r_i^2 $$

que es lo mismo que los residuos $r_i$ al cuadrado.

Y los mejores parámetros serán los que minimizen la suma de distancias cuadráticas:

$$ \sum_i \Big( y_i - (ax_i + b) \Big)^2 = \sum_i r_i^2 $$

Nota: por ahora, estoy asumiendo que los errores en $y_i$ son iguales.

En Python, pueden usar esta función:

def func(x, a, b):
    return a * x + b


parametros, covarianza = optimize.curve_fit(func, x, y)

que devuelve los parámetros $a$ y $b$:

parametros

array([ 1.82545875, 31.51416501])

y la matríz de covarianza:

covarianza

array([[ 0.00046089, -0.00668288],
       [-0.00668288,  0.13142996]])

que contiene la varianza (errores al cuadrado) de $a$ y $b$ en la diagonal:

np.sqrt(np.diag(covarianza))

array([0.02146831, 0.36253271])

fig, axes = plt.subplots(2, 1, figsize=(6, 8))

for ax in axes:
    ax.axhline(0, color="k")
    ax.axvline(0, color="k")
    ax.grid()
    ax.set(xlim=(-2, 32))

axes[0].set(ylim=(-5, 100))
axes[1].set(ylim=(-5, 5))

y_recta = func(x, *parametros)

axes[0].errorbar(x, y, yerr=1, fmt=".")
axes[0].plot(x, y_recta, color="r", linewidth=3)

axes[1].errorbar(x, y - y_recta, yerr=1, fmt=".")
axes[1].plot(x, y_recta - y_recta, "--r")

[<matplotlib.lines.Line2D at 0x7f1af2f73b80>]

Interpretación física, machine learning, y causalidad

El aprendizaje automático o machine learning está de moda en estos días (e Inteligencia Artificial y Big Data).

La respuesta está en los datos.

¿Pero está la respuesta ahí?

Regresión lineal (o ajuste lineal) es la base de machine learning.

Machine Learning y Estadistica están muy relacionados. Alguien (no tengo el link ahora) los diferenció por su objetivo:

• Estadística -> Estimar parámetros a partir de los datos: ¿cuáles son los $a$ y $b$ de mi modelo?

• Machine Learning -> Predecir a partir de los datos: ¿cuánto va a valer $y$ si observo $x$?

Entonces, ¿qué representan $a$ y $b$?

plt.scatter(x, y)
plt.plot(x, func(x, *parametros), color="r", linewidth=3)

plt.grid()
plt.axhline(0, color="k")
plt.axvline(0, color="k")
plt.xlim(-2, 32)
plt.ylim(-5, 100)

plt.xlabel("Temperatura [°C]")
plt.ylabel("Temperatura [°F]")

Text(0, 0.5, 'Temperatura [°F]')

Entonces, $a$ y $b$ son los coeficientes para cambiar de escalas de temperatura.

Pregunta tramposa:

¿Si yo pongo $x = 10 \text{ °C}$, voy a medir $y \approx 50 \text{ °F}$?

Es una pregunta causal. Para poder sacar esa conclusión, no solo hay que saber que representan las variables, sino como fue el experimento que las midió (en la jerga, se habla del data-generating process). ¿Cuál es el modelo causal que las conecta?.

Caso 1:

Había dos termómetros, uno en Celsius y el otro en Farenheit, y a lo largo de los días fuimos anotando los valores que mostraban cada uno.

Caso 2:

Había un termómetro en Celsius y un termostato en Farenheit, y a lo largo de los días fuimos anotando los valores que mostraban cada uno.

Periodo del péndulo

En la aproximación de ángulos pequeños:

$$ T^2 = \left( \frac{2\pi}{g} \right)^2 L $$

Ajustamos:

$$ y = a x + b $$

donde $ a = \left( \frac{2\pi}{g} \right)^2 $.

$b$ debería ser $0$. ¿Dió $b=0$? ¿A qué se puede deber?

$$\begin{align} T^2 &= \left( \frac{2\pi}{g} \right)^2 L + b \\\\&= a L + b \\\\&= a \Big(L - (-b/a) \Big) \\\\&= a \Big(L - L_0 \Big) \end{align}$$

Es decir, podemos interpretar la constante como un error al medir el largo del péndulo. Por ejemplo, no medir correctamente hasta el centro del masa. Hay que chequear que de un número razonable.

Extra:

Chequear que no se puede hacer esta interpretación si ajustabamos $T \propto \sqrt L$

Cuadrados mínimos

En general, podemos calcular la distancia cuadrática de los $y$ al modelo o función $f(x, p)$ que querramos:

$$ \Big( y_i - f(x_i, p) \Big)^2 = r_i^2 $$

donde $r_i$ son los residuos y $p$ los parámetros de nuestro modelo.

Y los mejores parámetros serán los que minimizen la suma de distancias cuadráticas:

$$ \sum_i \Big( y_i - f(x_i, p) \Big)^2 = \sum_i r_i^2 $$

En caso de que los $y_i$ no sean igual de confiables, podemos agregarle un peso $w_i$ a cada término:

$$ \sum_i w_i \; r_i^2 $$

Generalmente, se usa como peso la inversa de la varianza:

$$ w_i = \frac{1}{\sigma_i^2} $$

con lo que llegamos a:

$$ \chi^2 = \sum_i \left( \frac{y_i - f(x_i, p)}{\sigma_i} \right)^2 $$

Ejemplo "aburrido": función constante

Si $f(x, p) = p$, es decir, una constante que no depende de $x$, tenemos:

$$ \sum_i \Big( y_i - p \Big)^2 = \sum_i r_i^2 $$

datos = np.random.normal(loc=10, scale=1, size=30)


def func(y, p):
    return np.sum((y - p) ** 2)


t = np.linspace(7, 13, 100)
chi2 = np.array([func(datos, p) for p in t])


@ipywidgets.interact(p=(7, 13, 0.1))
def _(p=8):
    fig, ax = plt.subplots(1, 3, figsize=(12, 3))

    ax[0].set(ylim=(7, 13), title="Datos $y_i$")
    ax[1].set(ylim=(-3, 3), title="Residuos $r_i = y_i - p$")
    ax[2].set(xlabel="p", title="$\sum (y_i - p)^2$")

    ax[0].plot(datos, "o")
    ax[1].plot(datos - p, "o")
    ax[2].plot(t, chi2)

    ax[0].axhline(p, color="r")
    ax[1].axhline(0, color="k")
    ax[2].scatter(p, func(datos, p), color="r")

Podemos calcular que el mejor $p$, el que minimiza esa suma, es el promedio:

$$ p_{óptimo} = \bar{y} = \frac{1}{N} \sum_i y_i $$

A su vez, el promedio de los residuos para el $p$ óptimo es la varianza de los $y_i$:

$$ S^2 = \frac{1}{N} \sum_i \Big( y_i - \bar{y} \Big)^2 = \frac{1}{N} \sum_i r_i^2 $$

Es decir, a partir de los residuos podemos estimar los errores de los datos.

"Promedio que varia"

Cuadrados mínimos es una forma de calcular un promedio que varia en función de otra variable.

$$ \sum_i \Big( y_i - f(x_i, p) \Big)^2 = \sum_i \Big( y_i - \bar{y}(x_i, p) \Big)^2$$

x = np.arange(5)
for xi in x:
    y = np.random.normal(xi, scale=0.1, size=5)
    plt.plot(np.full_like(y, xi), y, "o")

plt.grid()

Ejemplo 2: funciones polinomiales

$$ f(x, p) = p_0 + p_1 \, x + p_2 \, x^2 + \ldots + p_N \, x^N = \sum_{i=0}^N p_i \, x^i $$

Supongamos que los datos provienen de una cuadrática con máximo (y raíz) en $x=2$:

$$ f(x) = 10 - (x - 2)^2 = -x^2 + 2 x - 6 $$

x = np.linspace(0, 4, 6)
y = 10 - (x - 2) ** 2
y = y + np.random.normal(scale=0.1, size=y.size)

plt.errorbar(x, y, yerr=0.1, fmt="--o")

<ErrorbarContainer object of 3 artists>

def func(x, a, b, c):
    return a * x ** 2 + b * x + c


popt, _ = optimize.curve_fit(func, x, y)
popt

array([-1.00603749,  4.04527684,  5.97677897])

Python incluye una función para hacer ajustes de polinomios de grado N:

np.polyfit(x, y, 2)

array([-0.96917382,  3.88894378,  6.11921903])

Cuadrática para refinar máximos

Puede ser útil ajustar una cuadrática para refinar la posición y el valor de un máximo local:

# Datos
plt.plot(x, y, "o", label="Datos")

# Ajuste
t = np.linspace(x[0], x[-1], 100)
plt.plot(t, func(t, *popt), label="Ajuste")

# Maximo
x_max = np.roots(np.polyder(popt))  # importa el order de los coeficientes
plt.plot(x_max, func(x_max, *popt), 'o', label=f"Maximo = {x_max[0]:.3f}")

plt.legend()

<matplotlib.legend.Legend at 0x7f1af2f16af0>

Sobreajuste

np.random.seed(1)

x = np.linspace(0, 6, 30)
y = 2 * x**3 + 3
y = y + np.random.normal(scale=1, size=y.size)

@ipywidgets.interact(grado=(1, y.size-1))
def _(grado=1):
    fig = plt.figure(figsize=(12, 4))
    left, right = fig.add_gridspec(1, 2)
    
    ax = left.subgridspec(2, 1).subplots(sharex=True)
    ax[0].plot(x, y, "o") # Datos

    # Ajuste
    coef = np.polyfit(x, y, grado)
    t = np.linspace(x[0], x[-1], 100)
    ax[0].plot(t, np.polyval(coef, t))

    # Residuos
    r = y - np.polyval(coef, x)
    ax[1].scatter(x, r)
    ax[1].axhline(0, color="k", ls="--")

    ax[0].set(title=f"$\sigma = {np.std(r):.2f}$")
#     ax[1].set(ylim=(-1.5, 1.5))
    
    ax_right = fig.add_subplot(right)
    ax_right.plot(x, y, "o")
    t = np.linspace(-5, 10, 100)
    ax_right.plot(t, np.polyval(coef, t))
    ax_right.set(title="Extrapolando")

    poly = " + ".join(f"{c:.2f} \, x^{i}" for i, c in enumerate(coef[::-1]))
    fig.suptitle(f"${poly}$", y=1.1)