Clase 13: Aprendizaje de Máquina (Supervisado)

Clase 13: Aprendizaje de Máquina (Supervisado)#

Recordemos: Aprendizaje a partir de datos previamente etiquetados $$x_i\to y_i.$$ Las etiquetas pueden ser categóricas o continuas (regresión).

Para el primer ejercicio vamos a hacer una regresión lineal para ajustar datos históricos de temperatura.

import pandas as pd
import requests
from io import StringIO

ESTACION = "80110099999"  # ID global de la estación (WMO) EOH en Medellin
AÑO_INICIO = 1965
AÑO_FIN = 2025

# Base URL del repositorio GSOD de la NOAA (formato CSV)
BASE_URL = "https://www.ncei.noaa.gov/data/global-summary-of-the-day/access/"

datos_completos = []

print("Iniciando descarga de datos desde la NOAA...")

for anio in range(AÑO_INICIO, AÑO_FIN + 1):
    url = f"{BASE_URL}{anio}/{ESTACION}.csv"

    try:
        # Realizamos la petición para el año específico
        response = requests.get(url, timeout=10)

        if response.status_code == 200:
            # Leemos el CSV directamente desde la respuesta en memoria
            df_anio = pd.read_csv(StringIO(response.text))
            datos_completos.append(df_anio)
            print(f"Año {anio}: Descargado con éxito.")
        elif response.status_code == 404:
            print(f"Año {anio}: No se encontraron datos (404).")
        else:
            print(f"Año {anio}: Error en el servidor (Código {response.status_code}).")

    except Exception as e:
        print(f"Año {anio}: Error en la conexión -> {e}")

# Concatenamos todos los años en un solo DataFrame
if datos_completos:
    df_total = pd.concat(datos_completos, ignore_index=True)

    # Convertir columna de fecha a formato datetime
    df_total['DATE'] = pd.to_datetime(df_total['DATE'])

    # --- PROCESAMIENTO Y CONVERSIÓN ---
    # La NOAA entrega la temperatura en Fahrenheit ('TEMP'). La convertimos a Celsius.
    df_total['TEMP_C'] = (df_total['TEMP'] - 32) * 5 / 9

    # Agrupamos por año y mes para obtener la media mensual que solicitaste
    df_total['AÑO_MES'] = df_total['DATE'].dt.to_period('M')

    # Extraemos el promedio mensual de la temperatura media, máxima y mínima
    # Nota: MAX y MIN también vienen originalmente en Fahrenheit.
    df_total['MAX_C'] = (df_total['MAX'] - 32) * 5 / 9
    df_total['MIN_C'] = (df_total['MIN'] - 32) * 5 / 9

    df_mensual = df_total.groupby('AÑO_MES').agg({
        'TEMP_C': 'mean',
        'MAX_C': 'mean',
        'MIN_C': 'mean'
    }).reset_index()

    # Guardar a un archivo CSV local
    df_mensual.to_csv("temperatura_mensual.csv", index=False)
    print(f"Archivo guardado como: 'temperatura_mensual.csv'")
    print(df_mensual.head())
else:
    print("\nNo se pudo descargar ningún dato.")

Iniciando descarga de datos desde la NOAA...

Año 1965: Descargado con éxito.

Año 1966: Descargado con éxito.

Año 1967: Descargado con éxito.

Año 1968: Descargado con éxito.

Año 1969: Descargado con éxito.

Año 1970: Descargado con éxito.

Año 1971: Descargado con éxito.

Año 1972: No se encontraron datos (404).

Año 1973: Descargado con éxito.

Año 1974: Descargado con éxito.

Año 1975: Descargado con éxito.

Año 1976: Descargado con éxito.

Año 1977: Descargado con éxito.

Año 1978: Descargado con éxito.

Año 1979: Descargado con éxito.

Año 1980: Descargado con éxito.

Año 1981: Descargado con éxito.

Año 1982: Descargado con éxito.

Año 1983: Descargado con éxito.

Año 1984: Descargado con éxito.

Año 1985: Descargado con éxito.

Año 1986: Descargado con éxito.

Año 1987: Descargado con éxito.

Año 1988: Descargado con éxito.

Año 1989: Descargado con éxito.

Año 1990: Descargado con éxito.

Año 1991: Descargado con éxito.

Año 1992: Descargado con éxito.

Año 1993: Descargado con éxito.

Año 1994: No se encontraron datos (404).

Año 1995: No se encontraron datos (404).

Año 1996: Descargado con éxito.

Año 1997: Descargado con éxito.

Año 1998: Descargado con éxito.

Año 1999: Descargado con éxito.

Año 2000: Descargado con éxito.

Año 2001: Descargado con éxito.

Año 2002: Descargado con éxito.

Año 2003: Descargado con éxito.

Año 2004: Descargado con éxito.

Año 2005: Descargado con éxito.

Año 2006: Descargado con éxito.

Año 2007: Descargado con éxito.

Año 2008: Descargado con éxito.

Año 2009: Descargado con éxito.

Año 2010: Descargado con éxito.

Año 2011: Descargado con éxito.

Año 2012: Descargado con éxito.

Año 2013: Descargado con éxito.

Año 2014: Descargado con éxito.

Año 2015: Descargado con éxito.

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[1], line 21
url = f"{BASE_URL}{anio}/{ESTACION}.csv"
try:
   # Realizamos la petición para el año específico
---> 21     response = requests.get(url, timeout=10)
   if response.status_code == 200:
       # Leemos el CSV directamente desde la respuesta en memoria
       df_anio = pd.read_csv(StringIO(response.text))

File /opt/homebrew/Caskroom/miniconda/base/envs/jupyter-book/lib/python3.10/site-packages/requests/api.py:73, in get(url, params, **kwargs)
def get(url, params=None, **kwargs):
   r"""Sends a GET request.

   :param url: URL for the new :class:`Request` object.
   (...)
   :rtype: requests.Response
   """
---> 73     return request("get", url, params=params, **kwargs)

File /opt/homebrew/Caskroom/miniconda/base/envs/jupyter-book/lib/python3.10/site-packages/requests/api.py:59, in request(method, url, **kwargs)
# By using the 'with' statement we are sure the session is closed, thus we
# avoid leaving sockets open which can trigger a ResourceWarning in some
# cases, and look like a memory leak in others.
with sessions.Session() as session:
---> 59     return session.request(method=method, url=url, **kwargs)

File /opt/homebrew/Caskroom/miniconda/base/envs/jupyter-book/lib/python3.10/site-packages/requests/sessions.py:592, in Session.request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json)
send_kwargs = {
   "timeout": timeout,
   "allow_redirects": allow_redirects,
}
send_kwargs.update(settings)
--> 592 resp = self.send(prep, **send_kwargs)
return resp

File /opt/homebrew/Caskroom/miniconda/base/envs/jupyter-book/lib/python3.10/site-packages/requests/sessions.py:706, in Session.send(self, request, **kwargs)
start = preferred_clock()
# Send the request
--> 706 r = adapter.send(request, **kwargs)
# Total elapsed time of the request (approximately)
elapsed = preferred_clock() - start

File /opt/homebrew/Caskroom/miniconda/base/envs/jupyter-book/lib/python3.10/site-packages/requests/adapters.py:645, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
   timeout = TimeoutSauce(connect=timeout, read=timeout)
try:
--> 645     resp = conn.urlopen(
       method=request.method,
       url=url,
       body=request.body,
       headers=request.headers,
       redirect=False,
       assert_same_host=False,
       preload_content=False,
       decode_content=False,
       retries=self.max_retries,
       timeout=timeout,
       chunked=chunked,
   )
except (ProtocolError, OSError) as err:
   raise ConnectionError(err, request=request)

File /opt/homebrew/Caskroom/miniconda/base/envs/jupyter-book/lib/python3.10/site-packages/urllib3/connectionpool.py:787, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
response_conn = conn if not release_conn else None
# Make the request on the HTTPConnection object
--> 787 response = self._make_request(
   conn,
   method,
   url,
   timeout=timeout_obj,
   body=body,
   headers=headers,
   chunked=chunked,
   retries=retries,
   response_conn=response_conn,
   preload_content=preload_content,
   decode_content=decode_content,
   **response_kw,
)
# Everything went great!
clean_exit = True

File /opt/homebrew/Caskroom/miniconda/base/envs/jupyter-book/lib/python3.10/site-packages/urllib3/connectionpool.py:464, in HTTPConnectionPool._make_request(self, conn, method, url, body, headers, retries, timeout, chunked, response_conn, preload_content, decode_content, enforce_content_length)
try:
   # Trigger any extra validation we need to do.
   try:
--> 464         self._validate_conn(conn)
   except (SocketTimeout, BaseSSLError) as e:
       self._raise_timeout(err=e, url=url, timeout_value=conn.timeout)

File /opt/homebrew/Caskroom/miniconda/base/envs/jupyter-book/lib/python3.10/site-packages/urllib3/connectionpool.py:1093, in HTTPSConnectionPool._validate_conn(self, conn)
# Force connect early to allow us to validate the connection.
if conn.is_closed:
-> 1093     conn.connect()
# TODO revise this, see https://github.com/urllib3/urllib3/issues/2791
if not conn.is_verified and not conn.proxy_is_verified:

File /opt/homebrew/Caskroom/miniconda/base/envs/jupyter-book/lib/python3.10/site-packages/urllib3/connection.py:796, in HTTPSConnection.connect(self)
   # Remove trailing '.' from fqdn hostnames to allow certificate validation
   server_hostname_rm_dot = server_hostname.rstrip(".")
--> 796     sock_and_verified = _ssl_wrap_socket_and_match_hostname(
       sock=sock,
       cert_reqs=self.cert_reqs,
       ssl_version=self.ssl_version,
       ssl_minimum_version=self.ssl_minimum_version,
       ssl_maximum_version=self.ssl_maximum_version,
       ca_certs=self.ca_certs,
       ca_cert_dir=self.ca_cert_dir,
       ca_cert_data=self.ca_cert_data,
       cert_file=self.cert_file,
       key_file=self.key_file,
       key_password=self.key_password,
       server_hostname=server_hostname_rm_dot,
       ssl_context=self.ssl_context,
       tls_in_tls=tls_in_tls,
       assert_hostname=self.assert_hostname,
       assert_fingerprint=self.assert_fingerprint,
   )
   self.sock = sock_and_verified.socket
# If an error occurs during connection/handshake we may need to release
# our lock so another connection can probe the origin.

File /opt/homebrew/Caskroom/miniconda/base/envs/jupyter-book/lib/python3.10/site-packages/urllib3/connection.py:975, in _ssl_wrap_socket_and_match_hostname(sock, cert_reqs, ssl_version, ssl_minimum_version, ssl_maximum_version, cert_file, key_file, key_password, ca_certs, ca_cert_dir, ca_cert_data, assert_hostname, assert_fingerprint, server_hostname, ssl_context, tls_in_tls)
   if is_ipaddress(normalized):
       server_hostname = normalized
--> 975 ssl_sock = ssl_wrap_socket(
   sock=sock,
   keyfile=key_file,
   certfile=cert_file,
   key_password=key_password,
   ca_certs=ca_certs,
   ca_cert_dir=ca_cert_dir,
   ca_cert_data=ca_cert_data,
   server_hostname=server_hostname,
   ssl_context=context,
   tls_in_tls=tls_in_tls,
)
try:
   if assert_fingerprint:

File /opt/homebrew/Caskroom/miniconda/base/envs/jupyter-book/lib/python3.10/site-packages/urllib3/util/ssl_.py:483, in ssl_wrap_socket(sock, keyfile, certfile, cert_reqs, ca_certs, server_hostname, ssl_version, ciphers, ssl_context, ca_cert_dir, key_password, ca_cert_data, tls_in_tls)
       context.load_cert_chain(certfile, keyfile, key_password)
context.set_alpn_protocols(ALPN_PROTOCOLS)
--> 483 ssl_sock = _ssl_wrap_socket_impl(sock, context, tls_in_tls, server_hostname)
return ssl_sock

File /opt/homebrew/Caskroom/miniconda/base/envs/jupyter-book/lib/python3.10/site-packages/urllib3/util/ssl_.py:527, in _ssl_wrap_socket_impl(sock, ssl_context, tls_in_tls, server_hostname)
   SSLTransport._validate_ssl_context_for_tls_in_tls(ssl_context)
   return SSLTransport(sock, ssl_context, server_hostname)
--> 527 return ssl_context.wrap_socket(sock, server_hostname=server_hostname)

File /opt/homebrew/Caskroom/miniconda/base/envs/jupyter-book/lib/python3.10/ssl.py:513, in SSLContext.wrap_socket(self, sock, server_side, do_handshake_on_connect, suppress_ragged_eofs, server_hostname, session)
def wrap_socket(self, sock, server_side=False,
               do_handshake_on_connect=True,
               suppress_ragged_eofs=True,
               server_hostname=None, session=None):
   # SSLSocket class handles server_hostname encoding before it calls
   # ctx._wrap_socket()
--> 513     return self.sslsocket_class._create(
       sock=sock,
       server_side=server_side,
       do_handshake_on_connect=do_handshake_on_connect,
       suppress_ragged_eofs=suppress_ragged_eofs,
       server_hostname=server_hostname,
       context=self,
       session=session
   )

File /opt/homebrew/Caskroom/miniconda/base/envs/jupyter-book/lib/python3.10/ssl.py:1104, in SSLSocket._create(cls, sock, server_side, do_handshake_on_connect, suppress_ragged_eofs, server_hostname, context, session)
       if timeout == 0.0:
           # non-blocking
           raise ValueError("do_handshake_on_connect should not be specified for non-blocking sockets")
-> 1104         self.do_handshake()
except (OSError, ValueError):
   self.close()

File /opt/homebrew/Caskroom/miniconda/base/envs/jupyter-book/lib/python3.10/ssl.py:1375, in SSLSocket.do_handshake(self, block)
   if timeout == 0.0 and block:
       self.settimeout(None)
-> 1375     self._sslobj.do_handshake()
finally:
   self.settimeout(timeout)

KeyboardInterrupt: 

import matplotlib.pyplot as plt
import seaborn as sns

df_mensual['AÑO_MES_DT'] = df_mensual['AÑO_MES'].dt.to_timestamp()

plt.figure(figsize=(15, 7))
sns.lineplot(x='AÑO_MES_DT', y='TEMP_C', data=df_mensual)
plt.xlabel('Date')
plt.ylabel('Average Temperature (°C)')
plt.grid(True)
plt.show()

import matplotlib.dates as mdates

# convert the datetime column to a numerical format
df_mensual['AÑO_MES_NUM'] = df_mensual['AÑO_MES_DT'].apply(mdates.date2num)

sns.lmplot(x='AÑO_MES_NUM', y='TEMP_C', data=df_mensual, ci=99, aspect=1.33, scatter_kws={'s': 10})

# make the x-axis readable
ax = plt.gca()
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
ax.xaxis.set_major_locator(mdates.YearLocator(5))
plt.gcf().autofmt_xdate()

plt.xlabel('Date')
plt.ylabel('Average Temperature (°C)')
plt.grid(True)
plt.show()

Solución exacta del problema de minimización - regresión lineal ideal#

Dado un conjunto de $m$ observaciones experimentales, queremos modelar la relación entre una variable dependiente e independiente a través de una línea recta (o un hiperplano en dimensiones mayores). El sistema sobredeterminado se plantea matricialmente como:

\[\begin{equation} A\theta \approx y \end{equation}\]

Para una regresión lineal simple con $n$ datos:

\[\begin{split}\begin{pmatrix} x_1 & 1 \\ x_2 & 1 \\ \vdots & \vdots \\ x_n & 1 \end{pmatrix} \begin{pmatrix} m_{\text{pendiente}} \\ b_{\text{intercepto}} \end{pmatrix} \approx \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix}\end{split}\]

$A$ (Matriz de diseño, dimensión $n \times 2$): Estructura que organiza nuestros datos conocidos. Cada fila es una observación. La primera columna contiene los datos independientes ($x_{\text{dato}}$) que ponderan a la pendiente, y la segunda columna es un vector de unos ($\mathbf{1}$) cuya función es mantener el intercepto constante para todas las ecuaciones.
$\theta$ (Vector de parámetros, dimensión $2 \times 1$): Es el vector de incógnitas ocultas que deseamos descubrir de forma óptima (la pendiente y el punto de corte).
$y$ (Vector de observaciones, dimensión $n \times 1$): Contiene los valores medidos de la variable dependiente.

Buscamos el estimador óptimo $\hat{\theta}$ que minimiza la norma euclidiana al cuadrado del residuo (la diferencia entre el modelo $A\theta$ y los datos reales $y$):

\[\begin{equation} \hat{\theta} = \arg\min_{\theta} \|A\theta - y\|_2^2 \end{equation}\]

La solución analítica exacta está dada por:

\[\hat{\theta} = \left(A^{T}A\right)^{-1}A^{T}y\]

Demostración#

Definimos la función de costo o funcional de error $J(\theta)$ como la norma $L_2$ al cuadrado del vector de residuos:

\[J(\theta) = \|A\theta - y\|_2^2\]

Utilizando las propiedades del producto punto y la transposición de matrices, expresamos la norma al cuadrado como el producto interno del vector de residuos consigo mismo:

\[J(\theta) = (A\theta - y)^T(A\theta - y)\]

Expandimos el producto distributivamente, manteniendo estrictamente el orden de las operaciones matriciales:

\[J(\theta) = (\theta^T A^T - y^T)(A\theta - y)\]

\[J(\theta) = \theta^T A^T A \theta - \theta^T A^T y - y^T A \theta + y^T y\]

Dado que el término $\theta^T A^T y$ resulta en un escalar (dimensión $1 \times 1$), es idéntico a su propia transpuesta. Al transponerlo obtenemos:

\[(\theta^T A^T y)^T = y^T (A^T \theta)^T = y^T A \theta\]

Como ambos términos centrales son equivalentes, los agrupamos en uno solo:

\[J(\theta) = \theta^T A^T A \theta - 2 y^T A \theta + y^T y\]

Para encontrar el mínimo global del funcional, calculamos el gradiente de $J(\theta)$ respecto al vector de parámetros $\theta$ ($\nabla_\theta$). Para ello, aplicamos las siguientes identidades del cálculo matricial:

$\nabla_\theta (\theta^T M \theta) = 2M\theta$ (válido debido a que la matriz kernel $M = A^T A$ es simétrica por definición: $(A^T A)^T = A^T A$).
$\nabla_\theta (v^T \theta) = v$, lo que transforma el término lineal en: $\nabla_\theta (2 (A^T y)^T \theta) = 2A^T y$.
El término $y^T y$ no depende de $\theta$, por lo que su derivada es cero.

Aplicando el gradiente término a término:

\[\nabla_\theta J(\theta) = 2A^T A \theta - 2A^T y\]

Igualamos el gradiente al vector nulo ($0$) para localizar el punto crítico que minimiza el error:

\[0 = 2A^T A \hat{\theta} - 2A^T y\]

Dividimos toda la expresión entre 2 y reordenamos los términos, llegando a las llamadas ecuaciones normales:

\[A^T A \hat{\theta} = A^T y\]

Finalmente, asumiendo que las columnas de la matriz de diseño $A$ son linealmente independientes, la matriz cuadrada $A^T A$ es de rango completo y, por lo tanto, invertible. Multiplicamos por la izquierda por su inversa $(A^T A)^{-1}$ para despejar el estimador óptimo:

\[\hat{\theta} = (A^T A)^{-1} A^T y\]

Mínimos Cuadrados Ponderados (WLS): Incorporando Errores#

Si cada observación $y_i$ tiene una incertidumbre conocida asociada $\sigma_i$ (desviación estándar del error), queremos minimizar la suma de los residuos normalizados por su propio error (el estadístico $\chi^2$):

\[\chi^2(\theta) = \sum_{i=1}^{m} \left( \frac{(A\theta)_i - y_i}{\sigma_i} \right)^2\]

Para operar esto matricialmente, definimos una matriz de pesos/ponderación diagonal $W$ de dimensión $n \times n$. Cada elemento de la diagonal es el recíproco de la varianza del error de ese dato:

\[\begin{split}W = \begin{pmatrix} 1/\sigma_1^2 & 0 & \cdots & 0 \\ 0 & 1/\sigma_2^2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1/\sigma_n^2 \end{pmatrix}\end{split}\]

Pre-multiplicamos todo nuestro sistema original por una matriz “raíz” $C$ (no hay que calcularla en el algoritmo, es sólo para la demostración), tal que $C^T C = W$. Como $W$ es diagonal, simplemente calculamos $C$ sacando la raíz cuadrada de cada término (es decir, el recíproco de la barra de error $1/\sigma_i$):

\[\begin{split}C = \begin{pmatrix} 1/\sigma_1 & 0 & \cdots & 0 \\ 0 & 1/\sigma_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 1/\sigma_n \end{pmatrix}\end{split}\]

Multiplicamos nuestro sistema original $A\theta \approx y$ por la izquierda por $C$, definiendo una nueva matriz de diseño modificada ($A_w$) y un nuevo vector de datos modificado ($y_w$):

\[\begin{split}A_w = CA = \begin{pmatrix} x_1/\sigma_1 & 1/\sigma_1 \\ x_2/\sigma_2 & 1/\sigma_2 \\ \vdots & \vdots \\ x_m/\sigma_n & 1/\sigma_n \end{pmatrix}, \quad y_w = Cy = \begin{pmatrix} y_1/\sigma_1 \\ y_2/\sigma_2 \\ \vdots \\ y_n/\sigma_mn \end{pmatrix}\end{split}\]

Como puedes ver, cada fila de la matriz de diseño y del vector $y$ se divide por su respectivo error $\sigma_i$.

Al aplicar la misma demostración analítica sobre el sistema modificado $A_w \theta \approx y_w$, la solución óptima resulta en:

\[\hat{\theta} = \left(A_w^T A_w\right)^{-1} A_w^T y_w\]

Sustituyendo $A_w = CA$ y $y_w = Cy$, y recordando que $C^T C = W$, la ecuación se reduce a la famosa fórmula de mínimos cuadrados ponderados:

\[\hat{\theta} = \left(A^T W A\right)^{-1} A^T W y\]

Ajuste Polinomial de Grado $n$ mediante Mínimos Cuadrados#

Aunque el modelo pase de una línea recta a una curva, el aparato matemático y la solución exacta de la demostración no cambian. El problema sigue siendo lineal respecto a los parámetros $\theta$.

La solución óptima sigue siendo: $$\hat{\theta} = (A^T A)^{-1} A^T y$$

El único cambio radica en la construcción de la Matriz de Diseño $A$ , la cual ahora incorpora las potencias de los datos independientes $x$.

Para ajustar un polinomio de grado $m$ del tipo $y = \theta_m x^m + \theta_{m-1} x^{m-1} + \dots + \theta_1 x + \theta_0$ a partir de $m$ datos, el sistema $A\theta \approx y$ se define como:

\[\begin{split}\begin{pmatrix} x_1^m & x_1^{m-1} & \cdots & x_1 & 1 \\ x_2^m & x_2^{m-1} & \cdots & x_2 & 1 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ x_n^m & x_n^{m-1} & \cdots & x_n & 1 \end{pmatrix} \begin{pmatrix} \theta_m \\ \theta_{m-1} \\ \vdots \\ \theta_1 \\ \theta_0 \end{pmatrix} \approx \begin{pmatrix} y_1 \\ y_2 \\ \vdots \\ y_n \end{pmatrix}\end{split}\]

import numpy as np
from scipy import stats

x = df_mensual['AÑO_MES_NUM'].values
y = df_mensual['TEMP_C'].values

A = np.vstack([x, np.ones(len(x))]).T


# Method 1: numpy.linalg.lstsq
# Returns: (coefficients, residuals, rank, singular_values)
coefficients_lstsq, residuals, rank, singular_values = np.linalg.lstsq(A, y, rcond=None)
slope_lstsq, intercept_lstsq = coefficients_lstsq
print(f"numpy.linalg.lstsq:  Slope = {slope_lstsq:.6f}, Intercept = {intercept_lstsq:.6f}")

# Method 2: Manual calculation using normal equation (beta = (A.T . A)^-1 . A.T . y)
# A is the design matrix (x with intercept), y is the target variable
b_manual = np.dot(np.dot(np.linalg.inv(np.dot(A.T, A)), A.T), y)
slope_manual, intercept_manual = b_manual
print(f"Manual (normal equation): Slope = {slope_manual:.6f}, Intercept = {intercept_manual:.6f}")

# Method 3: scipy.stats.linregress
# Directly returns (slope, intercept, r_value, p_value, std_err)
slope_linregress, intercept_linregress, r_value, p_value, std_err = stats.linregress(x, y)
print(f"scipy.stats.linregress: Slope = {slope_linregress:.6f}, Intercept = {intercept_linregress:.6f}")

%%timeit
np.linalg.lstsq(A, y, rcond=None)

%%timeit
np.dot(np.dot(np.linalg.inv(np.dot(A.T, A)), A.T), y)

%%timeit
stats.linregress(x, y)

# multiply by the number of days in a year (approx 365.25) and then by 10 for a decade.

# slope from numpy.linalg.lstsq
slope_degrees_per_day = slope_lstsq

days_in_a_year = 365.25
years_in_a_decade = 10

slope_degrees_per_decade = slope_degrees_per_day * days_in_a_year * years_in_a_decade

print(f"The temperature trend is {slope_degrees_per_decade:.4f} degrees Celsius per decade.")

standard_error_per_day = std_err

standard_error_per_decade = standard_error_per_day * days_in_a_year * years_in_a_decade

print(f"The standard error of the temperature trend is {standard_error_per_decade:.4f} degrees Celsius per decade.")

Ahora calculemos manualmente el intervalo de confianza (ver Apéndice). Para un ajuste lineal, la ecuación de las bandas de confianza describe una hipérbola.

slope_linregress, intercept_linregress, r_value, p_value, std_err = stats.linregress(x, y)

# calculate predicted y values (regression line)
y_pred = slope_linregress * x + intercept_linregress

# calculate confidence intervals
# residual standard error (sigma)
n = len(x)
residuals = y - y_pred
rmse = np.sqrt(np.sum(residuals**2) / (n - 2)) # RSE

# calculate sum of squared differences for x
x_mean = np.mean(x)
Sxx = np.sum((x - x_mean)**2)

# standard error of the mean response for each x_i
se_y_pred = rmse * np.sqrt(1/n + (x - x_mean)**2 / Sxx)

# critical t-value for 95% confidence interval
confidence_level = 0.95
degrees_freedom = n - 2
t_critical = stats.t.ppf((1 + confidence_level) / 2, degrees_freedom)

# calculate upper and lower bounds of the confidence interval
lower_bound = y_pred - t_critical * se_y_pred
upper_bound = y_pred + t_critical * se_y_pred

plt.figure(figsize=(15, 7))

plt.scatter(df_mensual['AÑO_MES_DT'], y, label='Monthly Average Temperature', s=10, alpha=0.6)

plt.plot(df_mensual['AÑO_MES_DT'], y_pred, color='red', label=f'Linear Regression (Slope: {slope_linregress*365.25*10:.4f} °C/decade)')

# confidence interval band
plt.fill_between(df_mensual['AÑO_MES_DT'], lower_bound, upper_bound, color='red', alpha=0.2, label='95% Confidence Interval')

ax = plt.gca()
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
ax.xaxis.set_major_locator(mdates.YearLocator(5))
plt.gcf().autofmt_xdate()

plt.xlabel('Date')
plt.ylabel('Average Temperature (°C)')
plt.grid(True)
plt.legend()
plt.show()

Cambiemos el rango de tiempo para ver lo que ocurre con la calidad del ajuste.

import matplotlib.dates as mdates
from scipy import stats

# filter df_mensual for the years 2005-2010
df_filtered = df_mensual[(df_mensual['AÑO_MES_DT'].dt.year >= 2005) & (df_mensual['AÑO_MES_DT'].dt.year <= 2008)].copy()

x_filtered = df_filtered['AÑO_MES_NUM'].values
y_filtered = df_filtered['TEMP_C'].values

slope_linregress_filtered, intercept_linregress_filtered, r_value_filtered, p_value_filtered, std_err_filtered = stats.linregress(x_filtered, y_filtered)

y_pred_filtered = slope_linregress_filtered * x_filtered + intercept_linregress_filtered

n_filtered = len(x_filtered)
residuals_filtered = y_filtered - y_pred_filtered
rmse_filtered = np.sqrt(np.sum(residuals_filtered**2) / (n_filtered - 2)) # Residual Standard Error

x_mean_filtered = np.mean(x_filtered)
Sxx_filtered = np.sum((x_filtered - x_mean_filtered)**2)

se_y_pred_filtered = rmse_filtered * np.sqrt(1/n_filtered + (x_filtered - x_mean_filtered)**2 / Sxx_filtered)

confidence_level = 0.95
degrees_freedom_filtered = n_filtered - 2
t_critical = stats.t.ppf((1 + confidence_level) / 2, degrees_freedom_filtered)

lower_bound_filtered = y_pred_filtered - t_critical * se_y_pred_filtered
upper_bound_filtered = y_pred_filtered + t_critical * se_y_pred_filtered

plt.figure(figsize=(15, 7))

plt.scatter(df_filtered['AÑO_MES_DT'], y_filtered, label='Monthly Average Temperature', s=10, alpha=0.6)

slope_degrees_per_decade_filtered = slope_linregress_filtered * 365.25 * 10
plt.plot(df_filtered['AÑO_MES_DT'], y_pred_filtered, color='red', label=f'Linear Regression (Slope: {slope_degrees_per_decade_filtered:.4f} °C/decade)')

plt.fill_between(df_filtered['AÑO_MES_DT'], lower_bound_filtered, upper_bound_filtered, color='red', alpha=0.2, label='95% Confidence Interval')

ax = plt.gca()
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m'))
ax.xaxis.set_major_locator(mdates.YearLocator(1)) # Show every year
plt.gcf().autofmt_xdate()

plt.xlabel('Date')
plt.ylabel('Average Temperature (°C)')
plt.grid(True)
plt.legend()
plt.show()

Concepto: Train / Validation / Test Sets#

Los datos suelen dividirse en subconjuntos con funciones distintas durante el entrenamiento y evaluación de modelos.

Training set#

Conjunto utilizado para ajustar el modelo.
El algoritmo aprende parámetros usando estos datos.
Ejemplo:
- pesos de una red neuronal,
- centroides en clustering,
- coeficientes de regresión.

Validation set#

Conjunto utilizado para evaluar configuraciones del modelo durante el desarrollo.
Se usa para:
- ajustar hiperparámetros,
- comparar modelos,
- detectar sobreajuste.
No debe utilizarse para el ajuste final de parámetros.
No siempre se usa.

Test set#

Conjunto reservado para evaluación final.
Simula datos completamente no vistos.
Se utiliza una sola vez al final para estimar desempeño real del modelo.
Si la regresión es continua, simplemente puede ser muestras de las variables independientes (p. ej. un np.linspace de un eje).

Conjunto	Propósito
Training	Aprender parámetros
Validation	Ajustar hiperparámetros y seleccionar modelos
Test	Evaluación final

Usualmente se toma la proporción train/test como 80/20.

Bootstrapping#

Técnica de remuestreo utilizada para estimar incertidumbre, estabilidad y error fuera de muestra (out-of-sample error).
Consiste en generar múltiples datasets re-muestreando aleatoriamente los datos originales con reemplazo.
Cada muestra bootstrap tiene el mismo tamaño que el dataset original, pero algunos puntos pueden repetirse y otros quedar excluidos.
Es como hacer muchas particiones train/test para estimar la robustez del modelo.

En validación de modelos#

Se genera una muestra bootstrap.
El modelo se ajusta usando esa muestra.
Se evalúa el desempeño sobre los datos no seleccionados (out-of-bag samples).
El proceso se repite muchas veces.

Out-of-bag (OOB)#

Aproximadamente un ( \sim 36% ) de los datos queda fuera de cada muestra bootstrap.
Esos datos actúan como conjunto de validación natural.
Permiten estimar error fuera de muestra sin separar explícitamente un test set fijo.

Método leave-one-out (jackknife)#

Se elimina un dato del conjunto.
El estimador se recalcula usando los datos restantes.
El proceso se repite para todos los datos.

Diferencia con bootstrap#

Método	Remuestreo
Bootstrap normal	Con reemplazo
Jackknife	Eliminando datos

Estudio de curvas de luz de transientes - MANTRA#

El proyecto MANTRA produjo curvas de luz de eventos transitorios astronómicos del Catalina Real-Time Transient Survey (CRTS), cubriendo 33000 sq deg del cielo desde 2007, usando datos de Mt. Lemmon Survey (MLS), Catalina Sky Survey (CSS), Siding Spring Survey (SSS). Cada uno de los 4869 eventos transitorios en MANTRA ha sido asociado a algún tipo de evento astronómico por un ser humano.

dfla=pd.read_csv("https://github.com/MachineLearningUniandes/MANTRA/raw/master/data/lightcurves/transient_labels.csv")
df=pd.read_csv("https://github.com/MachineLearningUniandes/MANTRA/raw/master/data/lightcurves/transient_lightcurves.csv")

Ejercicio: Ordene de mayor a menor los tipos de eventos transitorios en esta base de datos.

Inspeccionemos un par cualquiera de curvas de luz. Notemos que hay muchos datos faltantes, y que son distintos para distintos eventos. Adicionalmente, cada curva de luz tiene una cantidad distinta de datos. Esto significa que no podemos simplemente pasar todas las curvas por un solo algoritmo, sino que tenemos que hacer un pre-procesamiento intensivo antes.

from collections import Counter

ix=Counter(df.ID)
names=list(ix.keys())

i=320
name=names[i]
filt=df.ID==name
curve=(df[filt])
curve=curve.sort_values(by=['MJD'])
label=dfla.Classification[dfla.TransientID==float(name[6:])]
print(label,curve.shape)

plt.figure(figsize=(10, 6))
plt.plot(curve.MJD, curve.Mag, 'o', markersize=4)
plt.xlabel('Modified Julian Date (MJD)')
plt.ylabel('Magnitude')
plt.title(f'Light Curve for {label.iloc[0]} (ID: {name})')
plt.gca().invert_yaxis() # magnitudes
plt.show()

Gaussian Processes (GP)#

Modelo probabilístico no paramétrico para regresión y clasificación.
Un Gaussian Process es una colección de variables aleatorias indexadas por (x), tal que cualquier subconjunto finito tiene distribución gaussiana conjunta.
En lugar de aprender una función explícita, define una distribución sobre funciones.
Cualquier conjunto finito de puntos tiene una distribución conjunta gaussiana.
La inferencia se basa en una función de covarianza (kernel) que controla similitud entre puntos.
Si dos puntos, correspondientes a datos observados están cerca (según la escala del kernel), la covarianza es alta. Esto hace que la incertidumbre (varianza de la distribución posterior) para “interpolar” o “extrapolar” se reduce en el entorno de los puntos. Lejos de esos puntos, o si los puntos están lejos entre sí, la covarianza es baja y la incertidumbre en la predicción permanece alta.
El kernel incorpora supuestos sobre suavidad, periodicidad o escala de variación.
Tutorial interactivo: Enlace.

Visualización de ejemplo de GP#

El trabajo del Kernel: El rbf_kernel determina qué tan suaves son las líneas. Si disminuyes el length_scale, las líneas se volverán mucho más onduladas (con más curvas); si lo aumentas, se volverán mucho más suaves.
El Prior (Distribución a priori): Nota cómo el gráfico de la izquierda muestra curvas suaves, pero estas vagan sin rumbo fijo alrededor de una media de 0.
La Posterior (Distribución a posteriori): En el gráfico de la derecha, en el momento en que se le dan los puntos de entrenamiento rojos al GP (Proceso Gaussiano), todas las funciones aleatorias generadas se ven obligadas a pasar exactamente a través de esos puntos. Nota cómo la banda gris de incertidumbre se reduce a cero justo en los puntos de datos y se vuelve más ancha cuando se aleja de ellos.

import numpy as np
import matplotlib.pyplot as plt

# 1. Define the RBF Kernel (Your provided function)
def rbf_kernel(x1, x2, length_scale=1.0, sigma_f=1.0):
    """
    Squared Exponential (RBF) kernel
    """
    # Using np.subtract.outer or broadcasting to get pairwise squared distances
    sqdist = (x1[:, None] - x2[None, :])**2
    return sigma_f**2 * np.exp(-0.5 * sqdist / length_scale**2)

# 2. Setup the test inputs (the X axis where we want to plot our functions)
X_test = np.linspace(-5, 5, 100)

# Compute the covariance matrix for the test points (Prior Covariance)
K_ss = rbf_kernel(X_test, X_test, length_scale=1.0, sigma_f=1.0)

# The mean before seeing any data is just zeros
mu_prior = np.zeros(X_test.shape[0])

# --- PART 1: Sampling from the Prior ---
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
# Draw 3 random samples (functions) from the prior distribution
prior_samples = np.random.multivariate_normal(mu_prior, K_ss, 3)
for i, sample in enumerate(prior_samples):
    plt.plot(X_test, sample, label=f'Sample {i+1}')

plt.title("GP Prior (Before seeing data)")
plt.xlabel("X")
plt.ylabel("y")
plt.axhline(0, color='black', linestyle='--', alpha=0.3)
plt.legend()

# --- PART 2: Conditioning on Training Data (Posterior) ---
# Let's define 3 known data points (Training data)
X_train = np.array([-3.0, -1.0, 2.0])
y_train = np.array([1.5, -1.0, 0.5])

# Compute the rest of the required covariance matrices
K = rbf_kernel(X_train, X_train)
K_s = rbf_kernel(X_train, X_test)

# Apply the standard Gaussian Process conditioning formulas
# Mean: mu = K_s^T * K^-1 * y_train
# Covariance: Sigma = K_ss - K_s^T * K^-1 * K_s
K_inv = np.linalg.inv(K) # inv(K)

mu_posterior = K_s.T @ K_inv @ y_train
Sigma_posterior = K_ss - (K_s.T @ K_inv @ K_s)

plt.subplot(1, 2, 2)
# Draw 3 random samples from the posterior distribution
posterior_samples = np.random.multivariate_normal(mu_posterior, Sigma_posterior + 1e-6 * np.eye(100), 3)
for i, sample in enumerate(posterior_samples):
    plt.plot(X_test, sample, alpha=0.6)

# Plot the mean prediction and the actual data points
plt.plot(X_test, mu_posterior, 'k--', label='Mean Prediction', linewidth=2)
plt.scatter(X_train, y_train, color='red', zorder=5, s=50, label='Training Data')

# Plot the 95% confidence interval (2 standard deviations)
uncertainty = np.sqrt(np.diag(Sigma_posterior))
plt.fill_between(X_test, mu_posterior - 2*uncertainty, mu_posterior + 2*uncertainty, alpha=0.15, color='gray', label='95% Confidence')

plt.title("GP Posterior (After seeing data)")
plt.xlabel("X")
plt.legend()

plt.tight_layout()
plt.show()

Ahora vamos a usar una combinación de distintos tipos de kernels para ajustar nuestras curvas de luz.

Tendencia global (RBF): El componente RBF con escala de 1000 se encarga de modelar la tendencia macroscópica o de fondo de los datos (por ejemplo, el crecimiento lento a lo largo de meses).
Detalles locales (RationalQuadratic): El componente cuadrático racional se encarga de capturar los detalles, las fluctuaciones y los patrones más complejos o ruidosos a pequeña y mediana escala (por ejemplo, variaciones diarias).
Amplitud (ConstantKernel): Ajusta la escala vertical para que el rango de las predicciones coincida con el rango de los datos reales.

from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process import kernels as k

# Reshape MJD to a 2D array
X = curve.MJD.values.reshape(-1, 1)
y = curve.Mag.values

kernel = k.RationalQuadratic(length_scale_bounds=(1e-7, 1e5)) + k.ConstantKernel() + k.RBF(length_scale=1000)

# Initialize Gaussian Process Regressor
gpr = GaussianProcessRegressor(kernel=kernel, alpha=0.1, random_state=0)

# Fit the GPR model to the data
gpr.fit(X, y)

# Predict over a denser range of MJD values for a smooth curve
X_pred = np.linspace(X.min(), X.max(), 500).reshape(-1, 1)
y_mean_pred, y_std_pred = gpr.predict(X_pred, return_std=True)

# Plot the results
plt.figure(figsize=(12, 7))
plt.scatter(X, y, label='Original Data', s=15, alpha=0.7)
plt.plot(X_pred, y_mean_pred, color='red', label='GPR Mean Prediction')
plt.fill_between(X_pred.flatten(), y_mean_pred - 2 * y_std_pred, y_mean_pred + 2 * y_std_pred, color='red', alpha=0.2, label='95% Confidence Interval')

plt.xlabel('Modified Julian Date (MJD)')
plt.ylabel('Magnitude')
plt.title(f'Gaussian Process Regression for {label.iloc[0]} (ID: {name})')
plt.gca().invert_yaxis() # Invert y-axis for magnitudes
plt.grid(True)
plt.legend()
plt.show()

Apéndice - Propagación de Errores e Intervalos de Confianza#

Cuando las incertidumbres reales de la población ($\sigma$) son desconocidas, se estiman a posteriori utilizando la varianza de los residuos del ajuste. Al introducir esta estimación, la variable estandarizada adquiere una variabilidad adicional, lo que obliga a reemplazar la distribución Normal por la distribución $t$ de Student.

Supongamos que realizamos un ajuste polinomial o lineal utilizando la matriz de diseño $A$ (dimensión $m \times n$, donde $m$ es el número de observaciones y $n$ es el número de parámetros a estimar). Tras encontrar $\hat{\theta} = (A^T A)^{-1} A^T y$, calculamos el vector de residuos:

\[\hat{\epsilon} = y - A\hat{\theta}\]

La varianza residual del modelo ($\sigma^2$) se estima dividiendo la suma de los residuos al cuadrado entre los grados de libertad del sistema ($dof = m - n$):

\[s^2 = \frac{\hat{\epsilon}^T \hat{\epsilon}}{m - n} = \frac{\|y - A\hat{\theta}\|_2^2}{m - n}\]

La matriz de covarianza estimada de los parámetros ($\Sigma_\theta$) se calcula entonces escalando la matriz del sistema por esta varianza muestral:

\[\Sigma_\theta = s^2 (A^T A)^{-1}\]

Las incertidumbres estándar estimadas para cada parámetro individual corresponden a la raíz cuadrada de la diagonal principal: $s_{\theta_i} = \sqrt{(\Sigma_\theta)_{ii}}$.

Para construir el intervalo de confianza, analizamos la distribución estadística del estimador. Por el Teorema del Límite Central, sabemos que si los errores son normales, el estimador teórico estandarizado sigue una distribución normal:

\[Z = \frac{\hat{\theta}_i - \theta_i}{\sigma_{\theta_i}} \sim {N}(0, 1)\]

Por otro lado, la variable asociada a la estimación de la varianza muestral $s^2$ sigue una distribución Chi-cuadrado ($\chi^2$) dividida por sus grados de libertad:

\[V = \frac{(m-n)s^2}{\sigma^2} \sim \chi^2(m-n) \implies \frac{s^2}{\sigma^2} = \frac{V}{m-n}\]

Definimos una nueva variable aleatoria $T$ como el cociente entre la variable normal estándar $Z$ y la raíz cuadrada de la variable $\chi^2$ escalada. Matemáticamente, esta combinación define exactamente a la distribución $t$ de Student:

\[T = \frac{Z}{\sqrt{\frac{V}{m-n}}}\]

Sustituyendo las expresiones de $Z$ y del cociente de varianzas:

\[T = \frac{\frac{\hat{\theta}_i - \theta_i}{\sigma \sqrt{((A^TA)^{-1})_{ii}}}}{\sqrt{\frac{s^2}{\sigma^2}}} = \frac{\frac{\hat{\theta}_i - \theta_i}{\sigma \sqrt{((A^TA)^{-1})_{ii}}}}{\frac{s}{\sigma}} = \frac{\hat{\theta}_i - \theta_i}{s \sqrt{((A^TA)^{-1})_{ii}}}\]

Dado que el denominador es precisamente la incertidumbre estándar estimada ($s_{\theta_i}$), llegamos a la variable estandarizada final:

\[T = \frac{\hat{\theta}_i - \theta_i}{s_{\theta_i}} \sim t(m - n)\]

La variable $T$ sigue una distribución $t$ de Student con $\nu = m - n$ grados de libertad.

Definimos la probabilidad de que nuestra variable $T$ se encuentre confinada en una región simétrica delimitada por un valor crítico $t_{\text{crit}}$ para un nivel de confianza $1 - \alpha$:

\[P\left( -t_{\text{crit}} \le T \le t_{\text{crit}} \right) = 1 - \alpha\]

Sustituyendo la definición de $T$:

\[P\left( -t_{\text{crit}} \le \frac{\hat{\theta}_i - \theta_i}{s_{\theta_i}} \le t_{\text{crit}} \right) = 1 - \alpha\]

Multiplicamos toda la desigualdad por la incertidumbre estimada $s_{\theta_i}$:

\[P\left( -t_{\text{crit}} \cdot s_{\theta_i} \le \hat{\theta}_i - \theta_i \le t_{\text{crit}} \cdot s_{\theta_i} \right) = 1 - \alpha\]

Restamos el estimador calculado $\hat{\theta}_i$ en todos los miembros:

\[P\left( -\hat{\theta}_i - t_{\text{crit}} \cdot s_{\theta_i} \le - \theta_i \le -\hat{\theta}_i + t_{\text{crit}} \cdot s_{\theta_i} \right) = 1 - \alpha\]

Multiplicamos la desigualdad por $-1$ (invirtiendo el sentido de los operadores lógicos $\le$ a $\ge$):

\[P\left( \hat{\theta}_i + t_{\text{crit}} \cdot s_{\theta_i} \ge \theta_i \ge \hat{\theta}_i - t_{\text{crit}} \cdot s_{\theta_i} \right) = 1 - \alpha\]

Reordenando los límites de forma creciente, obtenemos el intervalo de confianza formal:

\[P\left( \hat{\theta}_i - t_{\text{crit}} \cdot s_{\theta_i} \le \theta_i \le \hat{\theta}_i + t_{\text{crit}} \cdot s_{\theta_i} \right) = 1 - \alpha\]

Esto demuestra que el parámetro real $\theta_i$ pertenece al intervalo $\hat{\theta}_i \pm t_{\text{crit}} \cdot s_{\theta_i}$ con una probabilidad exacta de $1-\alpha$.

Implementación en Python#

Debido a la simetría de la distribución $t$, el área acumulada bajo la curva desde $-\infty$ hasta el límite crítico de la derecha ($t_{\text{crit}}$) debe integrar el área de confianza central más la cola izquierda ($\alpha/2$):

\[\text{Área acumulada} = (1 - \alpha) + \frac{\alpha}{2} = \frac{2 - 2\alpha + \alpha}{2} = \frac{2 - \alpha}{2} = \frac{1 + (1 - \alpha)}{2}\]

Llamando $C = 1 - \alpha$ al nivel de confianza expresado en fracción (ej. 0.95), evaluamos la función cuantil o inversa de la distribución acumulada (PPF) de Student:

\[t_{\text{crit}} = F_t^{-1}\left(\frac{1 + C}{2}, \, \nu = m - n\right)\]

Lo que en código de Python se ejecuta de manera directa como:

from scipy import stats

# Grados de libertad = número de datos (m) - número de parámetros (n)
degrees_freedom = len(y_dato) - A.shape[1]

# Cálculo exacto del valor crítico t
t_critical = stats.t.ppf((1 + confidence_level) / 2, degrees_freedom)

Intervalo de confianza para un ajuste lineal#

Queremos predecir el valor de la recta en un punto arbitrario del eje independiente, al cual llamaremos $x_0$. La ecuación del modelo estimado es:

\[\hat{y}(x_0) = \hat{\theta}_1 x_0 + \hat{\theta}_0\]

Donde $\hat{\theta}_1$ es la pendiente y $\hat{\theta}_0$ es el intercepto calculados. Para construir el intervalo de confianza, necesitamos encontrar la varianza de esta predicción, $\text{Var}(\hat{y}(x_0))$.

Aplicando las propiedades del operador varianza para la suma de dos variables aleatorias dependientes (ya que la pendiente y el intercepto se calculan a partir del mismo conjunto de datos y están correlacionados):

\[\text{Var}(\hat{y}(x_0)) = \text{Var}(\hat{\theta}_1 x_0 + \hat{\theta}_0)\]

\[\text{Var}(\hat{y}(x_0)) = x_0^2 \cdot \text{Var}(\hat{\theta}_1) + \text{Var}(\hat{\theta}_0) + 2x_0 \cdot \text{Cov}(\hat{\theta}_1, \hat{\theta}_0)\]

Aquí queda explícita la naturaleza del problema: la varianza de la predicción depende directamente de $x_0^2$ y de $x_0$, lo que define matemáticamente una función cuadrática (una parábola en la varianza).

En un ajuste sin ponderación, la matriz de diseño $A$ para $m$ datos contiene la columna de la variable independiente y la columna de unos. La matriz de covarianza de los parámetros está dada por $\Sigma_\theta = s^2 (A^T A)^{-1}$. Al resolver esta inversión matricial analíticamente en términos de la media muestral $\bar{x} = \frac{\sum x_i}{m}$, obtenemos los siguientes componentes:

Varianza de la pendiente: $\text{Var}(\hat{\theta}_1) = \frac{s^2}{\sum (x_i - \bar{x})^2}$
Varianza del intercepto: $\text{Var}(\hat{\theta}_0) = s^2 \left( \frac{1}{m} + \frac{\bar{x}^2}{\sum (x_i - \bar{x})^2} \right)$
Covarianza entre ambos: $\text{Cov}(\hat{\theta}_1, \hat{\theta}_0) = -\frac{s^2 \bar{x}}{\sum (x_i - \bar{x})^2}$

Sustituyendo estos tres componentes en la ecuación de la varianza de la predicción y factorizando los términos comunes, la expresión se simplifica elegantemente a:

\[\text{Var}(\hat{y}(x_0)) = s^2 \left( \frac{1}{m} + \frac{(x_0 - \bar{x})^2}{\sum (x_i - \bar{x})^2} \right)\]

Dado que la incertidumbre estándar es la raíz cuadrada de la varianza, el intervalo de confianza formal para la recta en cualquier punto $x_0$ es:

\[\text{IC}(x_0) = \hat{y}(x_0) \pm t_{\text{crit}} \cdot s \sqrt{\frac{1}{m} + \frac{(x_0 - \bar{x})^2}{\sum (x_i - \bar{x})^2}}\]

Efecto de Hipérbola: El término cuadrático $(x_0 - \bar{x})^2$ bajo la raíz cuadrada genera que las bandas de confianza tengan una geometría hiperbólica.