Análisis Avanzado de Datos — 7. Métodos de kernel: estimación de densidades y regresión

estadistica
suavizacion
kernel
Estimación de densidades por kernel (KDE) y regresión por kernel (Nadaraya-Watson), el ancho de banda como control de flexibilidad, y su comparación con splines.
Author

Wilder Ramírez Delgado

Published

September 27, 2026

Análisis Avanzado de Datos — 7. Métodos de kernel: estimación de densidades y regresión

Open in Colab

👋 Sobre el autor

Wilder Ramírez Delgado convierte datos en decisiones. Científico de Datos, Arquitecto de IA y asesor tecnológico con más de 20 años de experiencia en analítica avanzada, machine learning, Big Data e IoT. Ingeniero Electrónico y Magíster en Analítica de Datos. CEO y fundador de Business Innovation Technology (BIT), asesora a entidades gubernamentales y empresas privadas en Colombia y es referente en analítica y normatividad para los sectores postal y logístico del país. Una serie pensada para ir de la teoría a la práctica, con ejemplos concretos en cada paso.

📌 Parte de la serie Análisis Avanzado de Datos.

Del knot a la ventana

En el artículo anterior de la serie construiste splines: elegías dónde poner los nodos (knots) y ajustabas polinomios por tramos que se conectan suavemente. Funciona muy bien, pero exige decidir una estructura. Ahora tomamos otro camino, todavía dentro de la suavización, con una idea más directa: estimar el comportamiento en cada punto a partir de los datos que están cerca, sin construir ninguna base polinómica ni colocar knots.

Esa idea sencilla —una ventana alrededor de cada punto y un peso que decae con la distancia— da lugar a dos herramientas: la estimación de densidades (¿cómo se distribuyen mis datos?) y la regresión por kernel (¿qué valor de \(y\) espero para cada \(x\)?). Y esa misma idea trae de vuelta, con otro nombre, el compromiso sesgo-varianza que ya conoces.

Objetivos de aprendizaje

Al terminar este notebook vas a entender cómo funcionan los métodos de suavización por kernel y cuándo conviene usarlos. En concreto, vas a ser capaz de:

  • Definir formalmente una densidad de probabilidad, un kernel y el estimador de densidad por kernel (KDE), y reconocer los tipos más usados (gaussiano, Epanechnikov, uniforme).
  • Explicar las limitaciones del histograma y estimar una densidad con el método de kernel (KDE), entendiendo su fórmula.
  • Aplicar la regresión por kernel (estimador de Nadaraya–Watson) para suavizar una relación no lineal.
  • Interpretar el papel del ancho de banda \(h\) como el parámetro que regula el compromiso sesgo-varianza, y seleccionarlo por validación.
  • Reconocer el sesgo de frontera y otras limitaciones prácticas, y comparar el suavizado por kernel con el spline y la recta.

Preparación del entorno

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import gaussian_kde
from sklearn.neighbors import KernelDensity
from sklearn.preprocessing import SplineTransformer
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import make_pipeline
from sklearn.model_selection import train_test_split, GridSearchCV, KFold
from sklearn.metrics import r2_score, mean_squared_error
from statsmodels.nonparametric.smoothers_lowess import lowess

SEMILLA = 2026
plt.rcParams['figure.figsize'] = (8, 4)
print("Librerías cargadas")
Librerías cargadas

Parte A — Estimación de densidades

Antes de suavizar una relación entre variables, empecemos por suavizar algo más básico: la forma en que se distribuye una sola variable. Es el terreno donde la idea del kernel se ve con más claridad.

Qué es una densidad (definición formal)

Cuando la variable es continua —aquí, la desviación en horas—, un valor exacto tiene probabilidad cero: \(P(X = -2.17) = 0\). Lo que sí tiene sentido es la probabilidad de un intervalo.

La función de densidad \(f\) es una función no negativa que convierte el área bajo la curva en probabilidad:

\[P(a \le X \le b) = \int_a^b f(x)\,dx, \qquad \int_{-\infty}^{\infty} f(x)\,dx = 1.\]

El número \(f(x)\) no es una probabilidad. Es cuánta probabilidad se concentra por unidad de \(x\) alrededor de ese punto (en este ejemplo, por hora). Por eso, en las figuras, los puntos negros van sobre el eje: marcan dónde cayeron los datos. La altura de la curva es densidad, no un \(y\) observado.

Un histograma de conteos responde «¿cuántos cayeron en este intervalo?». El mismo histograma con el eje en densidad —y, más adelante, el estimador por kernel— responde «¿qué tan concentrada está la probabilidad ahí?». El objetivo de esta parte es estimar \(f\) a partir de una muestra, sin suponer de antemano que sea una gaussiana.

El histograma engaña más de lo que crees

El histograma es la primera herramienta para ver una distribución, pero su forma depende de decisiones arbitrarias: el ancho de los intervalos (bins) y dónde empiezan. La misma variable puede contar historias distintas.

Generemos una variable con dos grupos. Piensa, por ejemplo, en la desviación (en horas) frente a la hora prometida de entrega para envíos de dos rutas mezclados: los valores negativos son entregas que llegaron antes; los positivos, entregas tardías. Vamos a dibujarla dos veces con el mismo ancho de bin, cambiando solo el punto donde empiezan los intervalos.

rng = np.random.default_rng(SEMILLA)
n = 600
comp = rng.random(n) < 0.5
datos = np.where(comp, rng.normal(-2, 0.6, n), rng.normal(2.5, 1.0, n))

ancho = 2.5
bordes_a = np.arange(-4.5, 8.0, ancho)      # intervalos que empiezan en -4.5
bordes_b = bordes_a + ancho / 4             # mismos intervalos, corridos 0.625

fig, ax = plt.subplots(1, 2, figsize=(11, 4), sharey=True)
conteos_a, _, _ = ax[0].hist(datos, bins=bordes_a, color='steelblue', edgecolor='white')
conteos_b, _, _ = ax[1].hist(datos, bins=bordes_b, color='steelblue', edgecolor='white')
ax[0].set_title(f'Ancho {ancho}, bordes desde {bordes_a[0]}')
ax[1].set_title(f'Ancho {ancho}, bordes desde {bordes_b[0]}')
for a in ax: a.set_xlabel('desviación frente a la hora prometida (horas)')
ax[0].set_ylabel('frecuencia')
plt.tight_layout(); plt.show()

print('Conteos, bordes desde -4.5  :', conteos_a.astype(int))
print('Conteos, bordes desde -3.875:', conteos_b.astype(int))
print('Media global:', round(datos.mean(), 2))

Conteos, bordes desde -4.5  : [153 149 210  88]
Conteos, bordes desde -3.875: [261  57 244  38]
Media global: 0.25

Cómo leer este resultado

Los dos gráficos usan exactamente los mismos datos y el mismo ancho de bin; solo se corrieron los bordes 0.625 horas. A la izquierda, los conteos son 153, 149, 210 y 88: se ve una meseta irregular, sin un valle claro. A la derecha, 261, 57, 244 y 38: aparecen dos grupos nítidos separados por un valle. La conclusión cambió por una decisión de dibujo que nada tiene que ver con los datos.

¿Por qué pasa? A la izquierda, un borde cae en −2, justo en el centro del primer grupo, y lo parte en dos barras de altura parecida. A la derecha, ese grupo cae casi entero dentro de un solo intervalo.

Fíjate además en la media global: 0.25, que cae en el valle entre los dos grupos, donde casi no hay envíos. Un resumen puede describir un punto que ningún dato ocupa. Para entender la variable necesitas su forma completa, y el histograma es demasiado frágil para dártela.

La misma trampa, en la avenida

En el artículo anterior el flujo de esa avenida se grafica contra la hora, y se ven dos picos: uno a las 8 y otro cerca de las 18. Ese gráfico conserva la hora. El histograma hace otra pregunta: si anotas solo la hora a la que pasa cada vehículo, ¿siguen viéndose las dos horas pico?

Abajo hay 12 000 vehículos repartidos con la misma curva de flujo: pasan más donde el flujo es más alto. Los dos paneles usan intervalos de 4 horas. Solo cambia dónde empiezan, a las 0 o a las 2. La línea roja es el flujo real, en vehículos por hora, para ver dónde cae cada corte. Las líneas punteadas marcan las 8:00 y las 18:12, el centro de cada pico.

def flujo_verdadero(hora):
    """Vehículos por minuto según la hora. La misma curva del artículo anterior."""
    ritmo_base = 24 + 7 * np.sin((hora - 5) * np.pi / 12)
    pico_manana = 80 * np.exp(-0.5 * ((hora - 8.0) / 0.95) ** 2)
    pico_tarde = 92 * np.exp(-0.5 * ((hora - 18.2) / 1.15) ** 2)
    return ritmo_base + pico_manana + pico_tarde


rng = np.random.default_rng(SEMILLA)
minutos = np.linspace(0, 24, 24 * 60, endpoint=False)
peso = np.clip(flujo_verdadero(minutos), 0, None)
peso = peso / peso.sum()
n_vehiculos = 12000
hora_paso = rng.choice(minutos, size=n_vehiculos, p=peso)

t = np.linspace(0, 24, 1000)
flujo = np.clip(flujo_verdadero(t), 0, None)
curva = n_vehiculos * flujo / np.trapezoid(flujo, t)

ancho = 4
inicios_a = list(range(0, 24, ancho))          # bordes en 0, 4, 8, 12, 16, 20
inicios_b = list(range(-2, 24, ancho))         # los mismos intervalos, corridos 2 horas

fig, ejes = plt.subplots(1, 2, figsize=(11, 4.2), sharey=True)
paneles = (
    (ejes[0], inicios_a, 'Ancho 4 h, bordes desde las 0'),
    (ejes[1], inicios_b, 'Ancho 4 h, bordes desde las 2'),
)
tasas = []
for eje, inicios, titulo in paneles:
    tasas_panel = []
    for inicio in inicios:
        lo = max(inicio, 0)
        hi = min(inicio + ancho, 24)
        if hi <= lo:
            continue
        conteo = np.sum((hora_paso >= lo) & (hora_paso < hi))
        tasa = conteo / (hi - lo)
        tasas_panel.append(round(tasa))
        eje.bar(
            lo, tasa, width=hi - lo, align='edge',
            color='steelblue', edgecolor='white', zorder=2,
        )
    eje.plot(t, curva, color='firebrick', lw=0.4, ls='-', zorder=3)
    eje.axvline(8.0, color='0.35', ls='--', lw=1, zorder=1)
    eje.axvline(18.2, color='0.35', ls='--', lw=1, zorder=1)
    eje.set_xlim(0, 24)
    eje.set_xticks(range(0, 25, 2))
    eje.set_title(titulo)
    eje.set_xlabel('hora del día')
    tasas.append(tasas_panel)

ejes[0].set_ylabel('vehículos por hora')
fig.legend(['flujo real'], loc='upper right', frameon=False, bbox_to_anchor=(0.99, 0.98))
fig.suptitle('Misma avenida, mismos vehículos: solo se corren los bordes 2 horas', y=1.03)
plt.tight_layout()
plt.show()

print('Bordes desde las 0:', tasas[0])
print('Bordes desde las 2:', tasas[1])

Bordes desde las 0: [206, 574, 642, 351, 958, 271]
Bordes desde las 2: [194, 261, 867, 361, 637, 673, 208]

Cómo leer este resultado

A la izquierda, con bordes en 0, 4, 8, 12, 16 y 20, las tasas son 206, 574, 642, 351, 958 y 271 vehículos por hora. El corte de las 8 parte la mañana en dos barras parecidas, 574 y 642. La tarde cae entera entre las 16 y las 20 y queda en 958, la barra más alta. Se lee un pico de tarde y una meseta de mañana, de 4 a 12.

A la derecha, con los bordes corridos 2 horas, los intervalos completos de 4 horas dan 261, 867, 361, 637 y 673. La mañana cae entera entre las 6 y las 10: 867, la barra más alta. El corte de las 18 parte la tarde en 637 y 673. Se lee lo contrario: un pico de mañana y una meseta de tarde, de 14 a 22.

Los extremos del panel derecho, 194 y 208, cubren solo dos horas (de 0 a 2 y de 22 a 24) porque el día se acaba. Su altura ya está en vehículos por hora, así que se comparan con las demás.

Los vehículos son los mismos. Cambia sobre qué pico cae el borde.

Qué es un kernel (definición formal)

Un kernel \(K\) es una función de ponderación: asigna a cada punto un peso según su distancia. En este artículo trabajaremos con kernels de suavización no negativos, normalizados y simétricos. Estas son las propiedades de los tres ejemplos siguientes; no son una definición universal para todos los métodos que usan la palabra kernel:

  • No negatividad: \(K(u) \ge 0\) para todo \(u\).
  • Integra uno: \(\int_{-\infty}^{\infty} K(u)\,du = 1\) (su área vale 1, como una densidad de probabilidad).
  • Simetría: \(K(-u) = K(u)\) (pesa igual a izquierda y derecha).

Hay muchos kernels válidos. Los tres más comunes:

  • Gaussiano: \(K(u) = \frac{1}{\sqrt{2\pi}} e^{-u^2/2}\). Suave, con soporte infinito (todo punto recibe algún peso).
  • Epanechnikov: \(K(u) = \frac{3}{4}(1-u^2)\) para \(|u|\le 1\). Tiene soporte finito y es óptimo en un sentido teórico preciso: entre los kernels no negativos y simétricos, y usando para cada uno su ancho de banda óptimo, minimiza el error integrado cuadrático medio asintótico (AMISE) en el marco habitual de estimación de densidades de segundo orden. Es un resultado asintótico bajo condiciones de regularidad, no una garantía de menor error en cualquier muestra.
  • Uniforme: \(K(u) = \frac{1}{2}\) para \(|u|\le 1\). Una caja: todos los vecinos dentro de la ventana pesan igual.
u = np.linspace(-3, 3, 400)
gauss = np.exp(-u**2 / 2) / np.sqrt(2 * np.pi)
epan = np.where(np.abs(u) <= 1, 0.75 * (1 - u**2), 0)
unif = np.where(np.abs(u) <= 1, 0.5, 0)

plt.plot(u, gauss, label='Gaussiano', linewidth=2)
plt.plot(u, epan, label='Epanechnikov', linewidth=2)
plt.plot(u, unif, label='Uniforme', linewidth=2)
plt.xlabel('u'); plt.ylabel('K(u)'); plt.legend()
plt.title('Tres kernels válidos'); plt.show()

Cómo leer este resultado

Los tres tienen la misma misión —dar más peso a lo cercano— pero con perfiles distintos. El gaussiano decae suavemente y nunca llega a cero: hasta los puntos lejanos aportan algo. El Epanechnikov y el uniforme tienen soporte finito: fuera de la ventana \(|u|\le 1\), el peso es exactamente cero. El uniforme es una caja (todos los vecinos pesan igual, como en un histograma móvil); el Epanechnikov es un compromiso suave con soporte acotado. Entre estos kernels habituales, la elección del ancho de banda suele influir más que la forma del kernel. La comparación teórica requiere optimizar el ancho para cada kernel y fijar un criterio de error; no significa que sean intercambiables con cualquier valor de h.

Qué es el KDE (definición formal)

La densidad \(f\) de la población no la observas. Observas una muestra \(x_1,\ldots,x_n\). El estimador de densidad por kernel (KDE, kernel density estimator) es una curva \(\hat f\) que estima esa \(f\) desconocida:

\[\hat{f}(x) = \frac{1}{n\,h} \sum_{i=1}^{n} K\!\left(\frac{x - x_i}{h}\right).\]

En palabras: colocas sobre cada dato una campana (el kernel \(K\)) de ancho \(h\), las sumas y divides por \(n\) para que el área siga siendo 1. \(\hat f(x)\) no es la densidad verdadera ni el valor de un dato: es una estimación, y cambia si cambian la muestra o el ancho de banda \(h\).

La diferencia con el histograma es que no partes el eje en intervalos fijos: cada observación lleva su propia campana centrada en ella.

Construyamos el KDE paso a paso

La muestra tiene 600 entregas. El KDE pone una campana sobre cada una, así que el estimador de verdad son 600 campanas promediadas. Si las dibujáramos todas a la vez, se fundirían y no verías el ladrillo. Por eso vamos en tres tomas, de menos a más:

  1. Un dato. Una sola campana, para ver la pieza básica: el kernel centrado en un punto.
  2. Ocho datos. Los primeros 8 de la muestra (datos[:8]), cada uno con su campana, sin sumarlas. Es una lupa: se distinguen los solapes y los huecos. No es el KDE todavía; es un recorte para que el ojo alcance.
  3. Los 600. Ahí sí se suman y se promedian todas las campanas. Esa curva ya es el KDE, en la misma escala que el histograma.

Para dibujar, no usamos los 600 datos como eje horizontal. Evaluamos cada campana en una malla de 500 valores de \(x\) entre \(-6\) y \(7\) horas (grid = np.linspace(-6, 7, 500)). Esas 500 columnas son puntos de dibujo, no observaciones: si pusieras 200 o 1000, la curva sería la misma idea, solo más tosca o más suave al trazarla.

Empezamos por la primera toma.

def campana(x, centro, h):
    """Kernel gaussiano centrado en 'centro' con ancho h (ya normalizado)."""
    return np.exp(-0.5 * ((x - centro) / h) ** 2) / (h * np.sqrt(2 * np.pi))

grid = np.linspace(-6, 7, 500)  # malla de dibujo: 500 valores de x, no 500 datos
punto = datos[0]
h = 0.3
altura = campana(punto, punto, h)

plt.plot(grid, campana(grid, punto, h), color='firebrick', linewidth=2)
plt.vlines(punto, 0, altura, colors='0.4', linestyles='--', linewidth=1)
plt.plot([punto], [0], 'ko', markersize=8, clip_on=False, zorder=5,
         label=f'un dato (x={punto:.2f})')
plt.xlim(punto - 2, punto + 2)
plt.ylim(-0.12, altura * 1.12)
plt.xlabel('desviación (horas)'); plt.ylabel('aporte'); plt.legend()
plt.title('Una campana (kernel) sobre un solo dato'); plt.show()

Cómo leer este resultado

Cada dato aporta una campana centrada en él. El punto negro vive sobre el eje —es el dato, no el valor de la densidad—. La línea punteada sube hasta el pico: ahí el aporte es máximo; a los lados, decae. Un solo punto genera un solo pico: la unidad básica con la que se construye toda la densidad.

Esta celda repite la campana, pero ya no sobre un solo dato: toma los primeros 8 de la muestra de entregas y dibuja una por cada uno.

  1. sub = datos[:8] se queda con esas ocho desviaciones, en horas.
  2. Cada punto negro es uno de esos datos, puesto sobre el eje en altura 0. El primer plot solo sirve para que la leyenda diga «datos» una vez; el bucle dibuja el resto.
  3. Para cada dato, campana(grid, xi, 0.3) es la misma gaussiana de antes, centrada en ese punto y con ancho \(h = 0.3\). Dividir por len(sub) (por 8) deja el aporte que ese dato tendría en el promedio: un octavo de la campana, no la campana entera. Por eso el eje se llama «aporte (escalado)».
  4. Las ocho curvas se dibujan por separado, sin sumarlas. Donde dos datos están cerca, las campanas se cruzan; donde hay un hueco, queda espacio vacío. En el siguiente paso ya no usamos este recorte: sumamos las campanas de los 600 datos.
# Los primeros 8 datos de la muestra de entregas.
sub = datos[:8]

# Un solo punto con etiqueta, para que la leyenda no repita «datos» ocho veces.
plt.plot([sub[0]], [0], 'ko', markersize=6, label='datos')
for xi in sub:
    # Campana centrada en xi, ancho 0.3, dividida por 8: el aporte de un dato al promedio.
    plt.plot(grid, campana(grid, xi, 0.3) / len(sub), color='steelblue', alpha=0.7)
    plt.plot([xi], [0], 'ko', markersize=6)
plt.xlabel('desviación (horas)'); plt.ylabel('aporte (escalado)'); plt.legend()
plt.title('Cada dato aporta su propia campana'); plt.show()

Cómo leer este resultado

Con varios puntos, cada uno pone su campana. Donde los datos se agolpan, las campanas se solapan y se apilan; donde hay huecos, apenas hay aporte. El siguiente paso es sumarlas todas.

Esta celda termina la construcción: suma las campanas de los 600 datos y las dibuja junto al histograma.

  1. Fija el ancho de cada campana en \(h = 0.3\) horas, el mismo de las figuras anteriores.
  2. Recorre cada dato y evalúa campana sobre la malla grid de 500 valores de \(x\) (la misma de las figuras anteriores: np.linspace(-6, 7, 500)). Sale una matriz de 600 filas × 500 columnas: una fila por entrega, una columna por cada punto donde se dibuja la curva. Las 500 no son datos; son el eje horizontal del gráfico.
  3. Promedia las filas. mean(axis=0) es \(\frac{1}{n}\) por la suma de las campanas. Como campana ya divide por \(h\), la curva integra uno: es una densidad. Esa curva es el KDE armado a mano.
  4. Dibuja el histograma con density=True, para que el eje sea densidad y no conteos, y superpone la curva. Así las dos estimaciones se leen en la misma escala.
# h: ancho de cada campana, en horas. El mismo de las figuras anteriores.
h = 0.3

# Una fila por cada uno de los 600 datos; 500 columnas = la malla de dibujo (grid), no 500 entregas.
aportes = np.array([campana(grid, xi, h) for xi in datos])

# Promedio de las 600 filas = (1/n) * suma. Eso es el KDE, y ya integra uno.
kde_manual = aportes.mean(axis=0)

# density=True pone el histograma en la misma escala que la curva (densidad, no conteos).
plt.hist(datos, bins=40, density=True, color='lightsteelblue',
         edgecolor='white', label='histograma')
plt.plot(grid, kde_manual, color='firebrick', linewidth=2, label='suma de campanas (KDE)')
plt.xlabel('desviación (horas)'); plt.ylabel('densidad'); plt.legend()
plt.title('El KDE es la suma (promediada) de todas las campanas'); plt.show()

Antes de seguir, pregúntate

Mirando la curva roja, ¿cuántos grupos (modas) muestra la variable? ¿Y esa forma dependió de dónde empezaban los intervalos, como pasaba con el histograma?

Cómo leer este resultado

Esa curva roja es el KDE: la suma de las 600 campanas, dividida por el número de datos para que integre uno. Cada campana está centrada en su propio dato, no en intervalos fijados de antemano, así que desaparece la arbitrariedad del punto de inicio que vimos en el histograma. Y la estructura real emerge sola: las dos modas que mezclamos. Acabas de construir un estimador de densidad desde cero. Lo que sí sigue por decidir es el ancho de cada campana, \(h\), que cumple el papel del ancho del bin; a eso volvemos enseguida.

En el ejemplo de las entregas, esos dos picos serían dos rutas: una que suele llegar unas dos horas antes de lo prometido, muy concentrada, y otra que llega unas dos horas y media tarde, más dispersa. La media global (0.25) no describe a ninguna de las dos: cae en el valle, donde casi no hay envíos. Mezclarlas en un solo número borra justo lo que te importa: que hay dos operaciones distintas.

La fórmula del KDE, ahora en código

La definición de arriba es lo que acabas de armar a mano:

\[\hat{f}(x) = \frac{1}{n\,h} \sum_{i=1}^{n} K\!\left(\frac{x - x_i}{h}\right)\]

donde \(K\) es el kernel, \(n\) el número de datos y \(h\) el ancho de banda. En la práctica no la programas a mano: scipy y scikit-learn la traen lista. Pero no hacen lo mismo con \(h\). scipy.gaussian_kde elige el ancho de banda por ti (regla de Scott por defecto, o Silverman si se lo pides con bw_method='silverman').

Las reglas de Scott y de Silverman son fórmulas cerradas para elegir \(h\) sin validar. En una dimensión las dos son del tipo \(h = c\,\hat\sigma\, n^{-1/5}\): miran el tamaño de muestra y la dispersión de los datos. El \(n^{-1/5}\) no es un número mágico: sale de minimizar el AMISE (asymptotic mean integrated squared error), un error teórico de toda la curva cuando \(n\) es grande —promedio, a lo largo de \(x\), de \((\hat f - f)^2\)—. Lo usan suponiendo que \(f\) se parece a una gaussiana (referencia normal). Scott, en scipy, usa el factor \(n^{-1/5}\) (\(c = 1\)). Silverman usa un factor distinto: en scipy es \((3n/4)^{-1/5}\); en la versión de los libros a veces entra también el rango intercuartílico, más resistente a valores extremos. Ninguna de las dos inspecciona si hay dos grupos o un pico raro: si la forma no es gaussiana, el \(h\) que entregan no es un óptimo para esa muestra, solo un atajo rápido. La fórmula del AMISE la desglosamos en la sección del compromiso sesgo-varianza.

En sklearn normalmente fijas \(h\) tú, y si le pides una regla automática hay un detalle que conviene conocer. Veamos esa diferencia.

kde_scipy = gaussian_kde(datos)  # bw_method='scott' por defecto
h_scott = float(np.sqrt(np.asarray(kde_scipy.covariance).squeeze()))

kde_sk_scott = KernelDensity(bandwidth='scott').fit(datos.reshape(-1, 1))
h_sk = float(kde_sk_scott.bandwidth_)
# score_samples devuelve log-densidad; exp la pasa a la misma escala que scipy.
dens_sk = np.exp(kde_sk_scott.score_samples(grid.reshape(-1, 1)))

print('scipy.gaussian_kde elige el ancho de banda solo (regla de Scott).')
print(f'  h de scipy:                  {h_scott:.3f}')
print(f'  h de sklearn con "scott":    {h_sk:.3f}')
print(f'  desviación estándar datos:   {datos.std(ddof=1):.3f}')

plt.hist(datos, bins=40, density=True, color='lightsteelblue',
         edgecolor='white', label='histograma')
plt.plot(grid, kde_scipy(grid), color='firebrick', linewidth=2,
         label=f'scipy (Scott, h={h_scott:.2f})')
plt.plot(grid, dens_sk, color='steelblue', linewidth=2,
         label=f'sklearn (Scott, h={h_sk:.2f})')
plt.xlabel('desviación (horas)'); plt.ylabel('densidad'); plt.legend()
plt.title('El mismo nombre Scott, dos anchos distintos'); plt.show()
scipy.gaussian_kde elige el ancho de banda solo (regla de Scott).
  h de scipy:                  0.660
  h de sklearn con "scott":    0.278
  desviación estándar datos:   2.373

Cómo leer este resultado

En el gráfico se ven las dos curvas con el mismo nombre «Scott». La roja (scipy, h = 0.660) es más ancha y aplana los dos grupos; la azul (sklearn, h = 0.278) se pega más a los datos y deja las dos modas más marcadas. No es que una librería «haga mejor Scott»: sklearn solo aplica el factor \(n^{-1/5}\) y no lo multiplica por la dispersión (desviación estándar de 2.373). \(0.278 \times 2.373 \approx 0.660\): el de scipy es el de sklearn por la escala de los datos. Solo equivalen si antes estandarizas la variable.

Fíjate en lo que no tuviste que decidir: el ancho de banda. scipy.gaussian_kde aplicó la regla de Scott que acabas de ver. Ese 0.660 es más ancho que el 0.3 que usamos a mano. Es esperable: estas reglas son de referencia normal, una aproximación asintótica basada en una forma gaussiana, no un óptimo garantizado para esta muestra. Con dos grupos, esa suposición falla y tienden a suavizar de más.

Por eso, con sklearn, lo habitual es fijar bandwidth=... explícitamente y elegirlo por validación, como haremos enseguida.

El ancho de banda y el compromiso sesgo-varianza

El KDE tiene un único parámetro crítico: el ancho de banda \(h\). Aquí el compromiso sesgo-varianza toma forma concreta:

  • Con \(h\) pequeño, cada campana es estrecha; la estimación se pega a cada dato y fluctúa mucho de una muestra a otra (varianza alta, sesgo bajo).
  • Con \(h\) grande, las campanas son anchas; la estimación es estable pero aplana los detalles reales (sesgo alto, varianza baja).

El error total se mide con el MISE (mean integrated squared error), \(\text{MISE}(h) = E\int (\hat f(x) - f(x))^2\,dx\). En una dimensión, para una densidad suficientemente suave y un kernel simétrico de segundo orden, con \(h \to 0\) y \(nh \to \infty\), su aproximación asintótica es

\[\text{AMISE}(h) \approx \underbrace{\frac{R(K)}{n\,h}}_{\text{varianza}} + \underbrace{\frac{h^4}{4}\,\mu_2(K)^2\,R(f'')}_{\text{sesgo}^2},\]

donde \(R(g) = \int g^2\) y \(\mu_2(K) = \int u^2 K(u)\,du\). El primer término baja cuando \(h\) crece; el segundo sube. El mínimo está en un punto intermedio, \(h^\ast \propto n^{-1/5}\), y de ahí sale el factor \(n^{-1/5}\) de la regla de Scott. Veámoslo con tres anchos.

anchos = [0.1, 0.3, 2.0]
fig, ax = plt.subplots(1, 3, figsize=(13, 4))
for a, hb in zip(ax, anchos):
    kde_h = KernelDensity(kernel='gaussian', bandwidth=hb).fit(datos.reshape(-1, 1))
    dens = np.exp(kde_h.score_samples(grid.reshape(-1, 1)))
    a.hist(datos, bins=40, density=True, color='lightsteelblue', edgecolor='white')
    a.plot(grid, dens, color='firebrick', linewidth=2)
    a.set_title(f'h = {hb}'); a.set_xlabel('desviación (horas)')
ax[0].set_ylabel('densidad')
plt.tight_layout(); plt.show()

Cómo leer este resultado

  • h = 0.1 (demasiado pequeño): curva ruidosa, con picos artificiales. Sobreajusta el ruido: varianza alta.
  • h = 0.3 (equilibrado): las dos modas se ven limpias, sin ruido espurio visible.
  • h = 2.0 (demasiado grande): las campanas anchas funden los dos grupos en una sola curva ancha, casi plana. Suaviza de más y borra la estructura: sesgo alto.

El principio es el mismo que con los knots del spline; solo cambia el parámetro que lo controla.

Elegir el ancho de banda por validación (densidad)

No queremos elegir \(h\) a simple vista. Como el KDE es una densidad, podemos medir qué tan bien predice datos nuevos mediante la log-verosimilitud en validación cruzada: buscamos el \(h\) que hace más probables los datos que el modelo no vio. GridSearchCV lo automatiza.

La verosimilitud es «qué tan creíbles se ven estos puntos si esta curva fuera la densidad verdadera»: se evalúa \(\hat f\) en cada dato y se multiplican esos valores. Como el producto de muchos números pequeños se va a cero, se toma el logaritmo y se suman: esa es la log-verosimilitud. Más alta = la curva coloca más densidad donde sí cayeron entregas nuevas. Más baja = las puso en el lugar equivocado. En validación cruzada se calcula con puntos que no se usaron para armar el KDE, para no premiar a un \(h\) que solo memoriza la muestra.

params = {'bandwidth': np.linspace(0.1, 2.0, 40)}
busqueda = GridSearchCV(KernelDensity(kernel='gaussian'), params, cv=5)
busqueda.fit(datos.reshape(-1, 1))
h_dens = busqueda.best_params_['bandwidth']
print(f'Ancho de banda óptimo por validación: {h_dens:.3f}')
Ancho de banda óptimo por validación: 0.295

Cómo leer este resultado

La validación elige un ancho de banda de 0.295, muy cerca del panel central (h = 0.3) que a simple vista se veía equilibrado, y menos de la mitad del 0.660 que proponía la regla de Scott. En este ejemplo, el criterio de validación favorece menos suavizado que la regla de Scott; no demuestra que ocurra así en toda distribución bimodal. Que el método basado en datos coincida con el juicio visual da confianza; además, a diferencia de la inspección visual, este procedimiento sigue un criterio explícito y reproducible. Ese es el \(h\) que usarías.

Parte B — Regresión por kernel

La misma idea de “mirar lo cercano” sirve para regresión. En vez de estimar una densidad, queremos estimar la función de regresión \(m(x) = E[Y \mid X = x]\), es decir, el valor esperado de \(y\) para cada \(x\). Lo haremos promediando los \(y\) de los puntos vecinos, con más peso para los más próximos.

La intuición: los pesos en un punto

Antes de la fórmula, veamos qué significa “promediar lo cercano”. Fijemos un punto \(x_0\) y miremos cuánto pesa cada dato según su distancia, usando el kernel gaussiano. Los cercanos deben pesar mucho; los lejanos, casi nada.

rng = np.random.default_rng(SEMILLA)
n = 200
x = np.sort(rng.uniform(0, 10, n))
m_true = np.sin(x) + 0.3 * x
y = m_true + rng.normal(0, 0.5, n)

x0 = 5.0
h = 0.5
pesos = np.exp(-0.5 * ((x - x0) / h) ** 2)

plt.scatter(x, y, c=pesos, cmap='Reds', s=30, edgecolor='gray', linewidth=0.3)
plt.axvline(x0, color='black', linestyle='--', label=rf'$x_0$ = {x0:g}')
plt.colorbar(label='peso del kernel')
plt.xlabel('x'); plt.ylabel('y'); plt.legend()
plt.title(r'Peso de cada dato al predecir en $x_0$ = 5'); plt.show()

Cómo leer este resultado

El color muestra el peso: los puntos junto a la línea \(x_0=5\) están en rojo intenso (pesan mucho), y el peso se desvanece rápido al alejarse. La predicción en \(x_0\) será un promedio de los \(y\), pero uno donde dominan los vecinos inmediatos. Cambiar \(x_0\) mueve la ventana; recorrer todos los \(x_0\) traza la curva suavizada.

El estimador de Nadaraya–Watson

Esa idea tiene fórmula: la predicción en \(x\) es el promedio de los \(y_i\) ponderado por el kernel según la distancia a \(x\):

\[\hat{m}(x) = \frac{\sum_{i=1}^{n} K\!\left(\frac{x - x_i}{h}\right) y_i}{\sum_{i=1}^{n} K\!\left(\frac{x - x_i}{h}\right)}\]

El numerador suma los \(y\) pesados; el denominador normaliza. No hay knots ni coeficientes: solo el ancho de banda \(h\). Lo implementamos a mano.

def _entradas_kernel(x_train, y_train, x_query, h):
    """Validar datos unidimensionales y un ancho de banda positivo."""
    x_train = np.asarray(x_train, dtype=float)
    y_train = np.asarray(y_train, dtype=float)
    x_query = np.atleast_1d(np.asarray(x_query, dtype=float))
    if x_train.ndim != 1 or y_train.ndim != 1 or x_query.ndim != 1:
        raise ValueError("Usa vectores unidimensionales; x_query también puede ser escalar.")
    if x_train.size == 0 or x_train.size != y_train.size:
        raise ValueError("x_train e y_train deben tener la misma longitud no nula.")
    if not all(np.isfinite(v).all() for v in (x_train, y_train, x_query)):
        raise ValueError("Los datos y las consultas deben ser finitos.")
    h = np.asarray(h, dtype=float)
    if h.ndim != 0 or not np.isfinite(h) or h <= 0:
        raise ValueError("h debe ser un escalar finito y positivo.")
    return x_train, y_train, x_query, float(h)


def _pesos_gaussianos(x_train, xq, h):
    # Restar el mayor log-peso evita que todos los pesos se redondeen a cero.
    # El factor común se cancela al normalizar y no cambia el estimador.
    with np.errstate(over="ignore", invalid="ignore"):
        log_w = -0.5 * ((x_train - xq) / h) ** 2
    if not np.isfinite(log_w).any():
        raise ValueError("La escala de distancias es extrema: reescala X o revisa h.")
    return np.exp(log_w - np.max(log_w))


def nadaraya_watson(x_train, y_train, x_query, h):
    x_train, y_train, x_query, h = _entradas_kernel(x_train, y_train, x_query, h)
    pred = np.empty(len(x_query))
    for i, xq in enumerate(x_query):
        w = _pesos_gaussianos(x_train, xq, h)
        pred[i] = np.sum((w / w.sum()) * y_train)
    return pred

grid_x = np.linspace(0, 10, 300)
ajuste = nadaraya_watson(x, y, grid_x, h=0.5)

plt.scatter(x, y, s=15, alpha=0.5, label='datos')
plt.plot(grid_x, np.sin(grid_x) + 0.3 * grid_x, 'k--', label='relación real')
plt.plot(grid_x, ajuste, color='firebrick', linewidth=2, label='kernel (h=0.5)')
plt.xlabel('x'); plt.ylabel('y'); plt.legend()
plt.title('Regresión por kernel (Nadaraya–Watson)'); plt.show()

Cómo leer este resultado

La curva roja sigue de cerca la relación real, \(m(x) = \sin(x) + 0.3x\) (línea negra), sin que le dijéramos que era un seno ni le ajustáramos un polinomio: solo promedió los \(y\) vecinos en cada punto. Eso es lo notable del suavizado por kernel: se adapta sin fijar una familia paramétrica global para la relación. El kernel y el ancho de banda sí introducen supuestos sobre la proximidad y la suavidad.

La implementación valida las entradas y estabiliza los pesos gaussianos para evitar divisiones 0/0 por redondeo numérico. Obtener una predicción finita lejos de los datos no vuelve confiable la extrapolación.

El ancho de banda controla la flexibilidad

Igual que en densidades, \(h\) decide cuánto suaviza. Cumple el papel que tenía el número de knots en el spline, pero en sentido inverso: más knots dan una curva más flexible; un \(h\) más grande da una curva más rígida.

anchos = [0.1, 0.5, 2.5]
plt.scatter(x, y, s=15, alpha=0.4, label='datos')
plt.plot(grid_x, np.sin(grid_x) + 0.3 * grid_x, 'k--', label='relación real')
colores = ['seagreen', 'firebrick', 'darkorange']
for hb, c in zip(anchos, colores):
    plt.plot(grid_x, nadaraya_watson(x, y, grid_x, hb), color=c, linewidth=2, label=f'h={hb}')
plt.xlabel('x'); plt.ylabel('y'); plt.legend()
plt.title('Efecto del ancho de banda'); plt.show()

Antes de seguir, pregúntate

¿Cuál de las tres curvas sigue el ruido punto por punto y cuál aplana las ondulaciones reales? ¿Con cuál te quedarías?

Cómo leer este resultado

  • h = 0.1 (verde): zigzaguea persiguiendo el ruido de cada dato. Sobreajuste.
  • h = 0.5 (rojo): sigue las ondulaciones reales sin copiar el ruido. Equilibrio.
  • h = 2.5 (naranja): aplana las ondas y casi vuelve a una recta. Subajuste.

Mismo compromiso: \(h\) pequeño = varianza alta; \(h\) grande = sesgo alto.

El sesgo de frontera (boundary bias)

Nadaraya–Watson tiene una debilidad conocida en los extremos del rango. En el interior, cada punto tiene vecinos a ambos lados y el promedio queda centrado. Pero en la frontera (el primer y el último tramo de \(x\)) solo hay vecinos de un lado. Si la relación tiene pendiente, esos vecinos están todos más arriba (o más abajo) y el promedio se sesga hacia el interior: la curva se “dobla”.

La solución clásica es cambiar el promedio local por una recta local. Nadaraya–Watson es una regresión polinómica local de grado 0: en cada punto ajusta una constante (el promedio ponderado). La regresión lineal local es de grado 1: en cada punto ajusta una recta por mínimos cuadrados ponderados con el mismo kernel, y usa su valor en ese punto como predicción.

Para que la comparación sea justa, usemos una relación lineal, \(y = 0.5x\) más ruido, y el mismo kernel gaussiano con el mismo \(h = 1\) en ambos métodos.

def lineal_local(x_train, y_train, x_query, h):
    x_train, y_train, x_query, h = _entradas_kernel(x_train, y_train, x_query, h)
    pred = np.empty(len(x_query))
    for i, xq in enumerate(x_query):
        w = _pesos_gaussianos(x_train, xq, h)  # mismo kernel gaussiano estabilizado
        # recta ponderada centrada en xq: el intercepto es la predicción en xq
        if np.unique(x_train[w > 0]).size < 2:
            raise ValueError("La recta local necesita dos valores distintos de X con peso positivo.")
        pendiente, intercepto = np.polyfit(x_train - xq, y_train, deg=1, w=np.sqrt(w))
        pred[i] = intercepto
    return pred

rng2 = np.random.default_rng(SEMILLA)
xb = np.sort(rng2.uniform(0, 10, 200))
mb = 0.5 * xb
yb = mb + rng2.normal(0, 0.6, 200)

gx = np.linspace(0, 10, 300)
fit_nw = nadaraya_watson(xb, yb, gx, h=1.0)
fit_ll = lineal_local(xb, yb, gx, h=1.0)

plt.scatter(xb, yb, s=12, alpha=0.35, label='datos')
plt.plot(gx, 0.5 * gx, 'k--', label='relación real (recta)')
plt.plot(gx, fit_nw, color='firebrick', linewidth=2,
         label='Nadaraya–Watson (grado 0)')
plt.plot(gx, fit_ll, color='seagreen', linewidth=2,
         label='lineal local (grado 1)')
plt.xlabel('x'); plt.ylabel('y'); plt.legend()
plt.title('Sesgo de frontera: promedio local vs recta local (mismo h)'); plt.show()

Cómo leer este resultado

Fíjate en los extremos. Nadaraya–Watson (rojo) se separa de la recta real cerca de \(x\approx 0\) y \(x\approx 10\): en el borde solo hay vecinos de un lado, el promedio se sesga hacia el interior y la curva se dobla. En el centro la sigue bien.

La lineal local (verde) usa el mismo kernel y el mismo \(h\); la única diferencia es que ajusta una recta en lugar de una constante. Esa pendiente extra le permite seguir la tendencia hasta el borde. Midamos cuánto mejora.

Cuánto se equivoca cada método en el borde

Como aquí conocemos la relación real, podemos medir el error de cada curva por zonas con la raíz del error cuadrático medio (RMSE, por sus siglas en inglés): los bordes (\(x < 1\) o \(x > 9\)) y el interior. Sumamos un tercer método, LOWESS de statsmodels: es una regresión lineal local con tres cambios prácticos (kernel tricúbico, una ventana definida por la fracción frac de vecinos más cercanos en lugar de un \(h\) fijo, e iteraciones que restan peso a los datos atípicos). Es la implementación de uso cotidiano de la idea de la recta local.

fit_low = lowess(yb, xb, frac=0.3, it=1, return_sorted=True)
fit_low_grid = np.interp(gx, fit_low[:, 0], fit_low[:, 1])

borde = (gx < 1) | (gx > 9)
m_grid = 0.5 * gx
ajustes = {'Nadaraya–Watson': fit_nw,
           'Lineal local':    fit_ll,
           'LOWESS':          fit_low_grid}

print('RMSE frente a la relación real')
print(f'{"método":<17}{"bordes":>8}{"interior":>10}')
for nombre, f in ajustes.items():
    rmse_borde = np.sqrt(np.mean((f[borde] - m_grid[borde]) ** 2))
    rmse_int = np.sqrt(np.mean((f[~borde] - m_grid[~borde]) ** 2))
    print(f'{nombre:<17}{rmse_borde:>8.3f}{rmse_int:>10.3f}')
RMSE frente a la relación real
método             bordes  interior
Nadaraya–Watson     0.283     0.091
Lineal local        0.020     0.038
LOWESS              0.058     0.065

Cómo leer este resultado

Los números confirman lo que se ve. En los bordes, Nadaraya–Watson tiene un RMSE de 0.283 frente a 0.091 en el interior: su error se triplica al llegar a los extremos. La lineal local, con el mismo kernel y el mismo \(h\), baja el error de borde a 0.020, y LOWESS a 0.058. Que en esta muestra la lineal local tenga incluso menos error en el borde que en el interior (0.038) es azar de la muestra, no una propiedad general.

Conviene separar sesgo y varianza. En este ejemplo la relación verdadera es exactamente lineal: si el ajuste está bien definido y el ruido tiene media condicional cero, la recta local reproduce esa tendencia sin sesgo condicional, incluso en el borde. El error observado en una sola muestra no es una medición directa del sesgo. Para una relación curva puede persistir sesgo de aproximación; además, disponer de información solo de un lado suele aumentar la incertidumbre. La recta local mejora el comportamiento en la frontera, pero no garantiza un menor error en toda muestra.

Elegir el ancho de banda por validación (regresión)

No queremos elegir \(h\) mirando el conjunto de prueba: dejaría de ser una medida honesta, porque habría participado en la decisión (una fuga de información). Reservamos los datos de prueba y elegimos \(h\) por validación cruzada (CV, por cross-validation) solo sobre el entrenamiento: el entrenamiento se divide en 5 partes (folds); cada una se usa por turnos para validar mientras las otras cuatro ajustan el modelo, y se promedia el error cuadrático medio (MSE). El conjunto de prueba solo se usa al final, para reportar.

X = x.reshape(-1, 1)
Xtr, Xte, ytr, yte, xtr, xte = train_test_split(
    X, y, x, test_size=0.3, random_state=SEMILLA
)

# K-Fold SOLO sobre el entrenamiento: el conjunto de prueba no entra a elegir h
kf = KFold(n_splits=5, shuffle=True, random_state=SEMILLA)
hs = np.linspace(0.05, 2.5, 50)
errores_cv = []
for h in hs:
    mse_folds = []
    for tr, va in kf.split(xtr):
        pred = nadaraya_watson(xtr[tr], ytr[tr], xtr[va], h)
        mse_folds.append(mean_squared_error(ytr[va], pred))
    errores_cv.append(np.mean(mse_folds))

h_opt = hs[int(np.argmin(errores_cv))]

plt.plot(hs, errores_cv, marker='o', color='steelblue')
plt.axvline(h_opt, color='firebrick', linestyle='--',
            label=f'h óptimo = {h_opt:.2f}')
plt.xlabel('ancho de banda h')
plt.ylabel('MSE de validación (K-Fold)')
plt.legend()
plt.title('Selección de h por CV (solo entrenamiento)')
plt.show()

print(f'h óptimo (CV en entrenamiento): {h_opt:.2f} | MSE de validación: {min(errores_cv):.3f}')

h óptimo (CV en entrenamiento): 0.30 | MSE de validación: 0.254

Cómo leer este resultado

Esta curva muestra el MSE de validación cruzada (promedio de 5 folds sobre el entrenamiento), no el error de prueba. El mínimo entre los candidatos está en h = 0.30, con un MSE de 0.254.

La varianza del ruido es \(\sigma^2 = 0.5^2 = 0.25\). Bajo este modelo generador y pérdida cuadrática, es el error irreducible en esperanza al predecir una nueva respuesta. No es una cota para cada MSE observado: por azar, una muestra puede producir un valor menor que 0.25, incluso al usar la función verdadera. La descomposición en ruido, sesgo al cuadrado y varianza describe promedios sobre nuevas observaciones y muestras de entrenamiento.

El valor 0.254 está cerca de esa referencia, pero no permite afirmar por sí solo que recuperamos casi toda la señal. Es el mínimo elegido entre varios candidatos y puede ser optimista por la propia selección; tampoco hemos cuantificado su incertidumbre. El conjunto de prueba reservado aporta la evaluación posterior.

A simple vista, \(h = 0.5\) parecía razonable. La CV favorece 0.30 dentro del rango explorado. El conjunto de prueba se usará para comparar recta, spline y kernel, sin volver a ajustar los parámetros a partir de ese resultado.

Comparación: recta, spline y kernel

Cerremos poniendo a competir los tres enfoques sobre los mismos datos de prueba, que no se usaron para elegir nada. El spline elige su número de knots por el mismo K-Fold sobre el entrenamiento; el kernel ya trajo su \(h\). Primero la gráfica; después el \(R^2\) en prueba.

candidatos_k = list(range(4, 16))
errores_k = []
for k in candidatos_k:
    mse_folds = []
    for tr, va in kf.split(Xtr):
        modelo_k = make_pipeline(
            SplineTransformer(n_knots=k, degree=3, include_bias=False),
            LinearRegression(),
        ).fit(Xtr[tr], ytr[tr])
        mse_folds.append(mean_squared_error(ytr[va], modelo_k.predict(Xtr[va])))
    errores_k.append(np.mean(mse_folds))
k_opt = candidatos_k[int(np.argmin(errores_k))]

recta = LinearRegression().fit(Xtr, ytr)
spline = make_pipeline(
    SplineTransformer(n_knots=k_opt, degree=3, include_bias=False),
    LinearRegression(),
).fit(Xtr, ytr)

r2_lineal = r2_score(yte, recta.predict(Xte))
r2_spline = r2_score(yte, spline.predict(Xte))
r2_kernel = r2_score(yte, nadaraya_watson(xtr, ytr, xte, h_opt))

plt.scatter(xte, yte, s=18, alpha=0.45, label='datos de prueba')
plt.plot(grid_x, recta.predict(grid_x.reshape(-1, 1)),
         color='gray', linewidth=2, label='recta')
plt.plot(grid_x, spline.predict(grid_x.reshape(-1, 1)),
         color='seagreen', linewidth=2, label=f'spline (k={k_opt})')
plt.plot(grid_x, nadaraya_watson(xtr, ytr, grid_x, h_opt),
         color='firebrick', linewidth=2, label=f'kernel (h={h_opt:.2f})')
plt.xlabel('x'); plt.ylabel('y'); plt.legend()
plt.title('Tres formas de ajustar la misma relación')
plt.show()

print(f'n_knots óptimo (CV en entrenamiento): {k_opt}')
print('R2 en prueba (no usado para elegir h ni knots)')
print(f'  Recta (lineal):     {r2_lineal:.3f}')
print(f'  Spline cúbico:      {r2_spline:.3f}')
print(f'  Kernel (NW, h opt): {r2_kernel:.3f}')

n_knots óptimo (CV en entrenamiento): 11
R2 en prueba (no usado para elegir h ni knots)
  Recta (lineal):     0.474
  Spline cúbico:      0.783
  Kernel (NW, h opt): 0.780

Cómo leer esta comparación

En el gráfico se ve de un vistazo lo que los números confirman. La recta gris atraviesa las ondas y se queda a media altura. El spline verde y el kernel rojo casi se superponen: ambos recuperan las ondulaciones sin copiar el ruido punto a punto.

  • Recta: 0.474. Deja más de la mitad de la variabilidad sin explicar.
  • Spline: 0.783 (con 11 knots elegidos por K-Fold en el entrenamiento) y Kernel: 0.780 (con h = 0.30, también por CV en el entrenamiento). Suben el \(R^2\) de 0.47 a cerca de 0.78: pasan de explicar menos de la mitad de la variabilidad a explicar casi cuatro quintas partes.

Importa el protocolo: \(h\) y el número de knots se eligieron sin mirar el conjunto de prueba, así que el \(R^2\) de prueba es una medida honesta. Pero es una medida con ruido: la prueba tiene solo 60 puntos, y una diferencia de 0.003 entre spline y kernel no permite declarar un ganador. La conclusión es que en esta partición muestran desempeños muy similares: no hay evidencia suficiente aquí para establecer una superioridad clara, pero tampoco se ha demostrado equivalencia estadística.

Que lleguen tan cerca no es casualidad: son dos caminos hacia el mismo objetivo, capturar una relación no lineal sin imponer de antemano una forma funcional. El spline usa una base de funciones y knots: una vez ajustado es un modelo lineal en esa base, predecir es barato y se integra de forma natural en modelos más grandes (como los GAM). El kernel promedia localmente: es más simple de entender y no exige decidir una estructura, pero necesita todos los datos de entrenamiento para cada predicción y sufre sesgo de frontera (la recta local lo reduce). Ninguno de los dos extrapola bien fuera del rango de los datos: fíjate en el tramo final del spline verde, a partir de \(x \approx 9.8\) (el último dato de entrenamiento), donde se vuelve plano porque SplineTransformer extrapola como una constante.

Errores comunes al trabajar con kernels

Síntoma Causa probable Solución
La densidad se ve llena de picos ancho de banda demasiado pequeño súbelo; elígelo por validación
Se pierde una moda o la asimetría ancho de banda demasiado grande bájalo; desconfía de las reglas automáticas con datos multimodales
La curva se dobla en los extremos sesgo de frontera usa regresión lineal local (LOWESS); interpreta los bordes con cautela
La curva se vuelve errática o queda indefinida donde hay pocos datos ventana sin vecinos suficientes; con soporte finito puede aparecer 0/0 y con un gaussiano los pesos pueden redondearse a cero limita el rango de predicción; estabiliza los pesos gaussianos y revisa \(h\) mediante validación
Resultados raros con varias variables en escalas distintas el kernel usa distancias estandariza las variables antes
bandwidth='scott' en sklearn da un \(h\) muy distinto al de scipy sklearn no reescala por la dispersión estandariza antes o fija bandwidth y valídalo
Muy lento con muchos datos el promedio recorre todos los puntos usa implementaciones optimizadas (árboles, FFT)

Lo que debes llevarte de este artículo

  • Los kernels de suavización usados aquí son funciones de ponderación no negativas, normalizadas y simétricas que favorecen la proximidad; hay varios (gaussiano, Epanechnikov, uniforme) y la elección importa poco frente al ancho de banda.
  • El KDE estima una densidad sumando una campana por dato; el estimador de Nadaraya–Watson estima una regresión promediando los \(y\) vecinos.
  • El ancho de banda \(h\) es el parámetro que controla el compromiso sesgo-varianza y se elige por validación, igual que el número de knots del spline (pero en sentido inverso: más \(h\), menos flexibilidad).
  • El promedio local sufre sesgo de frontera; la recta local (lineal local, LOWESS) lo reduce. Interpreta los extremos con cautela.
  • Spline y regresión por kernel abordan la estimación de una relación no lineal por caminos distintos. En este ejemplo rinden de forma similar; en otros datos el resultado debe comprobarse mediante evaluación.

¿Qué sigue?

Con esto cierras el tema de suavización: splines y kernels, dos formas de modelar relaciones no lineales sin imponer de antemano una forma funcional. La pregunta del próximo tema abre un frente nuevo:

¿Y si la variable que quieres predecir no es un número continuo, sino una categoría (sí/no) o un conteo?

Ahí entran los modelos lineales generalizados (GLM) —regresión logística, Poisson— y su extensión, los modelos aditivos (GAM), que combinan el GLM con suavizadores como los splines que ya conoces (y, en su versión clásica, también con suavizadores locales como los de este artículo). Todo lo que aprendiste sobre suavización reaparece ahí.

Referencias para profundizar

💬 ¿Te sirvió?

Deja en los comentarios una duda o un caso donde aplicarías esto — respondo todos. Sígueme para no perderte el próximo artículo de la serie y comparte con alguien que esté aprendiendo análisis de datos.

👉 El código completo está disponible para ejecutar directamente.