Más sobre ciencia de datos y estadística en: cienciadedatos.net

• Correlación lineal
• Ajuste y selección de distribuciones
• Kernel density estimation KDE
• Comparación de distribuciones con Python: test Kolmogorov–Smirnov
• Análisis de normalidad
• Análisis de homocedasticidad
• Test de permutación
• Bootstrapping
• T-test
• ANOVA.
• U-test

Introducción¶

Identificar el tipo de distribución que sigue una variable es un paso fundamental en prácticamente todos los estudios que implican datos, desde los contrastes de hipótesis hasta la creación de modelos de aprendizaje estadístico y machine learning.

Disponer de una función que describa aproximadamente los datos ofrece múltiples ventajas. Por ejemplo, permite calcular la densidad de probabilidad de una observación, estimar la probabilidad de que esta se encuentre dentro de un rango determinado, o simular nuevos valores sintéticos que sigan el mismo comportamiento que los originales.

En términos generales, ajustar una distribución consiste en encontrar la función matemática que mejor describa un conjunto de datos. De entre todas las funciones candidatas, el objetivo suele ser encontrar aquella que, con mayor probabilidad (verosimilitud), podría haber generado los datos observados.

Una de las aproximaciones más prácticas es utilizar distribuciones paramétricas. Estas son familias de distribuciones teóricas cuyo comportamiento está gobernado por un número fijo y finito de parámetros. Por ejemplo, la distribución normal queda totalmente determinada por su media ($\mu$) y su desviación típica ($\sigma$); ver más sobre este tipo de ajustes en Ajuste y selección de distribuciones con Python.

Sin embargo, cuando ninguna distribución paramétrica estándar se ajusta correctamente a los datos, es necesario recurrir a los métodos de ajuste no paramétricos. Su objetivo es generar una función de densidad que se adapte a la forma de los datos sin forzarlos a seguir una estructura rígida predefinida. Uno de los métodos más empleados es kernel density estimation (KDE).

En Python existen varias librerías que permiten ajustar distribuciones mediante KDE:

• SciPy: gaussian_kde.

• Statsmodels: KDEUnivariate y KDEMultivariate.

• Scikit-learn: KernelDensity.

A lo largo de este documento se describen los fundamentos teóricos del método KDE y ejemplos de cómo utilizar la implementación de Scikit-learn.

kernel density estimation (KDE)¶

En estadística, Kernel Density Estimation (KDE) es un método no paramétrico que permite estimar la función de densidad de probabilidad de una variable aleatoria continua a partir de un número finito de observaciones (muestra). Fue propuesto inicialmente por Fix y Hodges (1951) y desarrollado posteriormente por Rosenblatt (1956) y Parzen (1962).

Dado un valor $x$, la función aprendida por el kernel density estimator devuelve la densidad de la distribución en ese punto. Esta densidad, cuyo valor debe ser no negativo (rango $[0, +\infty)$), se interpreta como una medida relativa de verosimilitud (likelihood). Es importante notar que, en variables continuas, la densidad no es una probabilidad directa. Sin embargo, permite comparaciones: si la densidad en el punto A es mayor que en B, indica que es más probable encontrar observaciones en la vecindad de A que en la de B. Con frecuencia, para facilitar los cálculos numéricos y evitar problemas de desbordamiento, se utiliza el logaritmo de la densidad (log-likelihood); la interpretación se mantiene: cuanto mayor es su valor, mayor es la evidencia de que la observación es consistente con la distribución estimada.

Histograma¶

Una forma intuitiva de entender cómo funciona el KDE es partiendo del histograma, la herramienta más clásica para representar la distribución de datos. Su construcción sigue estos pasos:

1. Se identifica el rango de valores observados (mínimo y máximo).

2. Se divide ese rango en intervalos (bins), generalmente de igual tamaño.

3. Se cuenta el número de observaciones que caen dentro de cada intervalo.

4. Se representa una barra para cada intervalo, donde la altura es proporcional a la frecuencia (conteo).

5. Opcionalmente, se normaliza la altura de las barras para que el área total del histograma sume 1. En este caso, el eje $y$ representa la densidad de probabilidad empírica.

# Ejemplo histograma distribución bimodal
# ==============================================================================
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
plt.rcParams['lines.linewidth'] = 1.5
plt.rcParams['font.size'] = 8

muestra_1 = np.random.normal(loc=20, scale=5, size=200)
muestra_2 = np.random.normal(loc=40, scale=5, size=500)
datos     = np.hstack((muestra_1, muestra_2))

fig, ax = plt.subplots(figsize=(6, 3))
ax.hist(datos, bins=30, density=True, color="#3182bd", alpha=0.5)
ax.plot(datos, np.full_like(datos, -0.001), '|k', markeredgewidth=1)
ax.set_title('Histograma (normalizado)  \n 700 observaciones')
ax.set_xlabel('x')
ax.set_ylabel('densidad');

El principal inconveniente de los histogramas es que dependen en gran medida del número de intervalos (bins) que se utilicen, y de que la distribución generada es escalonada en lugar de continua. Este efecto se hace más evidente cuando hay pocas observaciones.

# Ejemplo histograma distribución bimodal pocas observaciones
# ==============================================================================
muestra_1 = np.random.normal(loc=20, scale=5, size=5)
muestra_2 = np.random.normal(loc=40, scale=5, size=5)
datos     = np.hstack((muestra_1, muestra_2))

fig, axs = plt.subplots(ncols=2, nrows=1, figsize=(8, 3))
axs[0].hist(datos, bins=30, density=True, color="#3182bd", alpha=0.5)
axs[0].plot(datos, np.full_like(datos, -0.001), '|k', markeredgewidth=1)
axs[0].set_title('Histograma (normalizado) \n 10 observaciones, 30 bins')
axs[0].set_xlabel('x')
axs[0].set_ylabel('densidad');

axs[1].hist(datos, bins=10, density=True, color="#3182bd", alpha=0.5)
axs[1].plot(datos, np.full_like(datos, -0.001), '|k', markeredgewidth=1)
axs[1].set_title('Histograma (normalizado) \n 10 observaciones, 10 bins')
axs[1].set_xlabel('x')
axs[1].set_ylabel('densidad');

Con solo 10 observaciones, incluso con diferente número de bins, es difícil que el histograma consiga una aproximación adecuada a la distribución real de la que proceden los datos.

Ajuste de un KDE¶

Idea intuitiva¶

kernel density estimation (KDE) expande la idea del histograma, "cada observación aumenta la densidad de probabilidad en la zona donde se encuentra", pero lo hace de forma que las contribuciones se agrupen creando una curva continua y suave (smooth). ¿Cómo lo consigue?

Supóngase que se dispone de un conjunto de observaciones $X$:

# Crear 4 observaciones
# ======================================================================================
X = np.array([2, 3, 4, 7])
fig, ax = plt.subplots(figsize=(5, 2.5))
ax.plot(X, np.full_like(X, 0.05), '|k', markeredgewidth=5)
ax.set_ylim(-0.05, 1)
ax.set_title('Posición observaciones')
ax.set_xlabel('x')
ax.set_ylabel('');

En lugar de colocar cada observación $x_i$ en un "cubo" rígido (como en el histograma), el método KDE sitúa una pequeña función matemática, llamada Kernel ($K$), centrada exactamente sobre cada uno de los puntos observados $x_i$. Uno de los kernel más utilizados es la distribución normal.

Siguiendo con el ejemplo, sobre cada observación se centra una distribución normal, con media igual al valor de la observación y desviación típica de 1 (más adelante se detalla la elección de este valor). De esta forma se consigue que cada observación contribuya justo en la posición que ocupa pero también, de forma gradual, en las regiones cercanas.

# Crear distribuciones normales centradas en cada observación
# ======================================================================================
from scipy import stats
fig, ax = plt.subplots(figsize=(5, 2.5))
ax.plot(X, np.full_like(X, 0), '|k', markeredgewidth=4)

Xgrid = np.linspace(min(X) - 1, max(X) + 1, num=500)
for x in X:
    densidad = stats.norm.pdf(Xgrid, loc=x, scale=1)
    ax.axvline(x=x, linestyle='--', color='black')
    ax.plot(Xgrid, densidad)
    
ax.set_title('Distribuciones normales (una por cada observación)')
ax.set_xlabel('x')
ax.set_ylabel('densidad');

Por último, si se suman las contribuciones individuales y se dividen por el total de curvas (observaciones), se consigue una curva final que describe la distribución de las observaciones.

# Agregar en cada punto la desidad de probabilidad de las 4 distribuciones
# ======================================================================================
fig, ax = plt.subplots(figsize=(5, 2.5))
ax.plot(X, np.full_like(X, 0), '|k', markeredgewidth=4)

Xgrid = np.linspace(min(X) - 1, max(X) + 1, num=500)
suma = np.full_like(Xgrid, 0)
for x in X:
    densidad = stats.norm.pdf(Xgrid, loc=x, scale=1)
    suma = suma + densidad
    ax.axvline(x=x, linestyle='--', color='black')
    ax.plot(Xgrid, densidad)
    
suma = suma / len(X)
ax.plot(Xgrid, suma, color='b')
ax.fill_between(Xgrid, suma, alpha=0.3, color='b')
ax.set_title('Suma de las distribuciones')
ax.set_xlabel('x')
ax.set_ylabel('densidad');

Esta idea tan sencilla pero a la vez tan potente es en la que se fundamenta el método kernel density estimation (KDE): aproximar una función de densidad como la suma de funciones (kernel) de cada observación.

Definición matemática¶

Dado un conjunto de datos $x = \{x_1, x_2, \ldots, x_n \}$, la función de densidad de probabilidad $f(x)$ puede aproximarse utilizando un kernel density estimation (KDE) tal que:

• $n$: es el número de datos (observaciones). Cada uno de ellos es el centro sobre el que se coloca un kernel.

• $h$: es el ancho de banda (bandwidth o smoothing parameter). Controla cuánto se expande la influencia de cada observación. Si se emplea como kernel una distribución normal, equivale a la desviación típica. Este es el valor más determinante a la hora de ajustar un KDE, puesto que condiciona el nivel de sobreajuste.

• $K$: es el Kernel, una función simétrica y positiva que integra a 1, que define la forma y la distribución de la influencia (peso) que se asocian a cada observación. En los ejemplos anteriores se ha utilizado como kernel la distribución normal.

Selección del ancho de banda¶

El ancho de banda es crucial a la hora de estimar una función densidad mediante el método KDE. Si su valor es muy bajo, se genera overfitting y la función resultante estará demasiado influenciada por el "ruido" de los datos. Si su valor es muy elevado, la función resultante no será capaz de aprender la distribución subyacente.

fig, axs = plt.subplots(ncols=3, nrows=1, figsize=(9, 3))
X = np.array([2, 3, 4, 7])
Xgrid = np.linspace(min(X) - 1, max(X) + 1, num=500)

axs[0].plot(X, np.full_like(X, 0), '|k', markeredgewidth=4)
suma = np.full_like(Xgrid, 0)
for x in X:
    densidad = stats.norm.pdf(Xgrid, loc=x, scale=10)
    suma = suma + densidad
suma = suma / len(X)
axs[0].plot(Xgrid, suma, color='b')
axs[0].fill_between(Xgrid, suma, alpha=0.3, color='b')
axs[0].set_title('Kernel: gausiano\n Ancho de banda: 10')
axs[0].set_xlabel('x')
axs[0].set_ylabel('densidad');

axs[1].plot(X, np.full_like(X, 0), '|k', markeredgewidth=4)
suma = np.full_like(Xgrid, 0)
for x in X:
    densidad = stats.norm.pdf(Xgrid, loc=x, scale=1)
    suma = suma + densidad 
suma = suma / len(X)
axs[1].plot(Xgrid, suma, color='b')
axs[1].fill_between(Xgrid, suma, alpha=0.3, color='b');
axs[1].set_title('Kernel: gausiano \n Ancho de banda: 1')
axs[1].set_xlabel('x')

axs[1].plot(X, np.full_like(X, 0), '|k', markeredgewidth=4)
suma = np.full_like(Xgrid, 0)
for x in X:
    densidad = stats.norm.pdf(Xgrid, loc=x, scale=0.1)
    suma = suma + densidad 
suma = suma / len(X)
axs[2].plot(Xgrid, suma, color='b')
axs[2].fill_between(Xgrid, suma, alpha=0.3, color='b');
axs[2].set_title('Kernel: gausiano \n Ancho de banda: 0.1')
axs[2].set_xlabel('x');

Aunque no existe una fórmula cerrada para conocer de antemano el valor óptimo del ancho de banda, existen estrategias consolidadas para identificar valores adecuados.

Reglas Empíricas (Heurísticas)

Son fórmulas rápidas que asumen cierta suavidad en los datos subyacentes. Las más comunes para kernels gaussianos son:

• Scott’s rule: $h \approx 1.06 \cdot \hat{\sigma} n^{-1/5}$

• Silverman’s rule: $h = 0.9 \cdot \min ( \hat{\sigma}, IQR/1.35 ) n^{-1/5}$

Si bien son computacionalmente inmediatas, estos métodos tienden a suavizar en exceso (oversmoothing) cuando la distribución real de los datos es multimodal o muy diferente a una normal.

Validación cruzada

Es el método más riguroso, ya que no asume ninguna forma predefinida para los datos. Se basa en probar distintos valores de $h$ y medir qué tan bien generalizan en datos no vistos.

Este método requiere mayor coste computacional, pero es aplicable a cualquier distribución. En Scikit-learn, esto se implementa habitualmente mediante GridSearchCV, utilizando la log-verosimilitud (log-likelihood) como métrica de evaluación. Cuando el número de observaciones es muy bajo, es altamente recomendable utilizar la estrategia Leave-One-Out Cross-Validation (LOOCV) para maximizar el uso de la información disponible en cada iteración.

Tipo de kernel¶

El kernel es la función que determina cómo se distribuye la influencia de cada observación, por lo tanto, puede tener un impacto notable en la estimación de la función de densidad resultante. Aunque en la gran mayoría de casos se emplea un kernel gaussiano (distribución normal), existen otras posibilidades. Las implementadas en scikit learn son:

• Gaussian: asigna los pesos siguiendo la distribución normal con una desviación estándar equivalente al ancho de banda.

• Epanechnikov: las observaciones que están a una distancia entre 0 y h tienen un peso entre $\frac{3}{4}$ y 0 con disminución cuadrática. Toda observación fuera de este rango tiene peso 0.

• Tophat: Asigna el mismo peso a todas las observaciones que estén dentro del ancho de banda.

• Exponential: el peso decae de forma exponencial.

• Linear: el peso decae de forma lineal dentro del ancho de banda. Más allá de este el peso es 0.

• Cosine: el peso dentro del ancho de banda es proporcional al coseno.

# Gráficos de los kernels disponibles (fuente: documentación scikit learn)
# ----------------------------------------------------------------------
from sklearn.neighbors import KernelDensity
X_plot = np.linspace(-6, 6, 1000)[:, None]
X_src = np.zeros((1, 1))

fig, ax = plt.subplots(2, 3, sharex=True, sharey=True, figsize=(8,4))
fig.subplots_adjust(left=0.05, right=0.95, hspace=0.05, wspace=0.05)

def format_func(x, loc):
    if x == 0:
        return '0'
    elif x == 1:
        return 'h'
    elif x == -1:
        return '-h'
    else:
        return '%ih' % x

for i, kernel in enumerate(['gaussian', 'tophat', 'epanechnikov',
                            'exponential', 'linear', 'cosine']):
    axi = ax.ravel()[i]
    log_dens = KernelDensity(kernel=kernel).fit(X_src).score_samples(X_plot)
    axi.fill(X_plot[:, 0], np.exp(log_dens), '-k', fc='#AAAAFF')
    axi.text(-2.6, 0.95, kernel)

    axi.xaxis.set_major_formatter(plt.FuncFormatter(format_func))
    axi.xaxis.set_major_locator(plt.MultipleLocator(1))
    axi.yaxis.set_major_locator(plt.NullLocator())

    axi.set_ylim(0, 1.05)
    axi.set_xlim(-2.9, 2.9)

ax[0, 1].set_title('Kernels disponibles Scikit-learn');

Ejemplo univariante¶

En el siguiente ejemplo se utiliza el método kernel density estimation (KDE) para tratar de ajustar la función de densidad de unos datos simulados a partir de una distribución normal bimodal (solapamiento de dos distribuciones normales).

Librerías¶

Las librerías utilizadas en este ejemplo son:

# Tratamiento de datos
# ==============================================================================
import pandas as pd
import numpy as np

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
plt.rcParams['lines.linewidth'] = 1.5
plt.rcParams['font.size'] = 8

# Ajuste de distribuciones
# ==============================================================================
from scipy import stats
from sklearn.neighbors import KernelDensity
from sklearn.model_selection import GridSearchCV
from sklearn.model_selection import KFold

# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('once')

Datos¶

Para esta demostración, se emplean datos simulados a partir de una distribución mixta formada por dos distribuciones normales:

• Distribución 1: $\mathcal{N}(1,\,0.5)$, peso 0.25

• Distribución 2: $\mathcal{N}(-1,\,0.5)$, peso 0.75

Los pesos indican la proporción en que cada distribución individual contribuye a la distribución final.

# Datos
# ==============================================================================
n = 1000
np.random.seed(123)
muestra_1 = np.random.normal(loc=1, scale=0.5, size=int(n * 0.75))
muestra_2 = np.random.normal(loc=-1, scale=0.5, size=int(n * 0.25))
datos = np.hstack((muestra_1, muestra_2))

Como los datos han sido simulados, se puede calcular la verdadera curva de densidad.

# Superposición de histograma con la verdadera distribución
# ==============================================================================
X_grid = np.linspace(-3, 4, 1000)
densidad = (stats.norm.pdf(loc=1, scale=0.5, x=X_grid)*0.75
              + stats.norm.pdf(loc=-1, scale=0.5, x=X_grid)*0.25)

fig, ax = plt.subplots(figsize=(6, 3))
ax.hist(datos, bins=30, density=True, color="#3182bd", alpha=0.5)
ax.plot(datos, np.full_like(datos, -0.001), '|k', markeredgewidth=1)
ax.plot(X_grid, densidad, color = 'blue', label='Distribución real')
ax.set_title('Distribución observada vs real')
ax.set_xlabel('x')
ax.set_ylabel('densidad')
ax.legend();

Modelo KDE¶

La clase KernelDensity de Scikit learn implementa el algoritmo para estimar la función de densidad de distribuciones univariantes y multivariantes.

En primer lugar, se entrena un modelo utilizando un kernel lineal y un ancho de banda (bandwidth) de 1.

# Ajuste del modelo KDE
# ==============================================================================
modelo_kde = KernelDensity(kernel='linear', bandwidth=1)
modelo_kde.fit(X=datos.reshape(-1, 1))

KernelDensity(bandwidth=1, kernel='linear')

Una vez entrenado el modelo, con el método score_samples() se puede predecir la densidad estimada para cualquier nuevo valor. El valor devuelto está en escala logarítmica.

# Predicción de densidad
# ==============================================================================
new_X = np.array([1])
log_density_pred = modelo_kde.score_samples(X=new_X.reshape(-1, 1))

# Se aplica el exponente para deshacer el logaritmo
density_pred = np.exp(log_density_pred)
density_pred

array([0.46014262])

Si se calcula la predicción para todo el rango de valores observado, se puede representar la función de densidad obtenida.

# Gráfico distribución de densidad estimada (predicción) vs real
# ==============================================================================
X_grid = np.linspace(-3, 4, 1000)
densidad_real = (
    stats.norm.pdf(loc=1, scale=0.5, x=X_grid) * 0.75 
    + stats.norm.pdf(loc=-1, scale=0.5, x=X_grid) * 0.25
)

log_densidad_pred = modelo_kde.score_samples(X_grid.reshape((-1, 1)))
#Se aplica el exponente para deshacer el logaritmo
densidad_pred = np.exp(log_densidad_pred)

fig, ax = plt.subplots(figsize=(6, 3))
ax.hist(datos, bins=30, density=True, color="#3182bd", alpha=0.5)
ax.plot(datos, np.full_like(datos, -0.001), '|k', markeredgewidth=1)
ax.plot(X_grid, densidad_real, color = 'blue', label='Distribución real')
ax.plot(X_grid, densidad_pred, color = 'red', label='KDE')
ax.set_title('Comparación de distribuciones')
ax.set_xlabel('x')
ax.set_ylabel('densidad')
ax.legend();

Selección de kernel y bandwidth¶

El tipo de kernel y bandwidth son hiperparámetros cuyo valor óptimo no puede conocerse de antemano. Se emplea validación cruzada para tratar de identificar la combinación que consigue mejores resultados.

Cuando se combina GridSearchCV con KernelDensity el score utilizado es el total log-likelihood. Cuanto más próximo a cero, mejor se ajusta el KDE a los datos.

# Validación cruzada para identificar kernel y bandwidth
# ==============================================================================
param_grid = {
    'kernel': ['gaussian', 'epanechnikov', 'exponential', 'linear'],
    'bandwidth': np.logspace(-2, 3, 20)
}

grid = GridSearchCV(
        estimator  = KernelDensity(),
        param_grid = param_grid,
        n_jobs     = -1,
        cv         = 10, 
        verbose    = 0
      )
_ = grid.fit(X = datos.reshape((-1, 1)))

resultados = (
    pd.DataFrame(grid.cv_results_)
    .loc[:, ['param_bandwidth', 'param_kernel', 'mean_test_score']]
    .sort_values(by='mean_test_score', ascending=False)
)
resultados.head(10)

# Mejores hiperparámetros por validación cruzada
# ==============================================================================
print("----------------------------------------")
print("Mejores hiperparámetros encontrados (cv)")
print("----------------------------------------")
print(grid.best_params_, ":", grid.best_score_)

modelo_kde_final = grid.best_estimator_

----------------------------------------
Mejores hiperparámetros encontrados (cv)
----------------------------------------
{'bandwidth': np.float64(0.37926901907322497), 'kernel': 'epanechnikov'} : -136.081077189232

# Gráficos distribución de densidad modelo final
# ==============================================================================
X_grid = np.linspace(-3, 4, 1000)
densidad_real = (
    stats.norm.pdf(loc=1, scale=0.5, x=X_grid) * 0.75
    + stats.norm.pdf(loc=-1, scale=0.5, x=X_grid) * 0.25
)

log_densidad_pred = modelo_kde_final.score_samples(X_grid.reshape((-1, 1)))
# Se aplica el exponente para deshacer el logaritmo
densidad_pred = np.exp(log_densidad_pred)

fig, ax = plt.subplots(figsize=(6, 3))
ax.hist(datos, bins=30, density=True, color="#3182bd", alpha=0.5)
ax.plot(datos, np.full_like(datos, -0.001), '|k', markeredgewidth=1)
ax.plot(X_grid, densidad_real, color = 'blue', label='Distribución verdadera')
ax.plot(X_grid, densidad_pred, color = 'red', label='Kernel: epanechnikov \n bw:0.34')
ax.set_title('Comparación de distribuciones')
ax.set_xlabel('x')
ax.set_ylabel('densidad')
ax.legend();

La función de densidad obtenida empleando un kernel de epanechnikov con un ancho de banda de 0.38 consigue aproximarse bien a la distribución real.

Cabe recordar que esta comparación solo es posible porque los datos han sido simulados. En la práctica, no se suele conocer la verdadera distribución, de lo contrario no se necesitaría aplicar una aproximación KDE.

Ejemplo multivariante¶

Todo lo explicado sobre el ajuste mediante kernel density estimation (KDE) aplica a distribuciones de más de una dimensión (multivariante). Sin embargo, es importante tener en cuenta que, debido al problema conocido como curse of dimensionality, los resultados empeoran exponencialmente a medida que se aumenta el número de variables.

Librerías¶

Las librerías utilizadas en este documento son:

# Tratamiento de datos
# ==============================================================================
import pandas as pd
import numpy as np

# Gráficos
# ==============================================================================
import matplotlib.pyplot as plt
import seaborn as sns
plt.style.use('fivethirtyeight')
plt.rcParams['lines.linewidth'] = 1.5
plt.rcParams['font.size'] = 8

# Ajuste de distribuciones
# ==============================================================================
from scipy import stats
from sklearn.neighbors import KernelDensity
from sklearn.model_selection import GridSearchCV, LeaveOneOut
from sklearn.preprocessing import StandardScaler

# Configuración warnings
# ==============================================================================
import warnings
warnings.filterwarnings('once')

Datos¶

El set de datos geyser contiene información sobre las erupciones del géiser Old Faithful localizado en el parque nacional de Yellowstone, Wyoming. En concreto, recoge información sobre la duración de 272 erupciones, así como el intervalo de tiempo transcurrido entre ellas.

# Lectura de datos
# ==============================================================================
datos = sns.load_dataset("geyser")
datos = datos[['duration', 'waiting']]
print(datos.info())
datos.head(3)

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 272 entries, 0 to 271
Data columns (total 2 columns):
 #   Column    Non-Null Count  Dtype  
---  ------    --------------  -----  
 0   duration  272 non-null    float64
 1   waiting   272 non-null    int64  
dtypes: float64(1), int64(1)
memory usage: 4.4 KB
None

# Histograma de cada variable
# ==============================================================================
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(8, 3))

axs[0].hist(datos['duration'], bins=30, density=True, color="#3182bd", alpha=0.5)
axs[0].set_title('Distribución de duración')
axs[0].set_xlabel('duración (duration)')
axs[0].set_ylabel('densidad')

axs[1].hist(datos['waiting'], bins=30, density=True, color="#3182bd", alpha=0.5)
axs[1].set_title('Distribución intervalo entre erupciones')
axs[1].set_xlabel('intervalo (waiting)')
axs[1].set_ylabel('densidad');

# Gráfico
# ==============================================================================
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(4, 4))
ax.scatter(datos['duration'], datos['waiting'], color="#3182bd", alpha=0.5)
ax.set_title('Distribución erupciones geiser')
ax.set_xlabel('duración (duration)')
ax.set_ylabel('intervalo (waiting)');

Escalado de variables¶

Cuanto se trabaja con distribuciones multivariantes, es fundamental escalar las variables antes de ajustar un KDE. En la implementación de Scikit-learn, KernelDensity utiliza un único ancho de banda para todas las dimensiones junto con una métrica de distancia euclidiana por defecto. Como resultado, si las variables tienen diferentes escalas numéricas o unidades, las dimensiones con rangos más grandes dominarán el cálculo de distancias y, por lo tanto, la forma de la densidad estimada. Esto hace que el KDE refleje diferencias en la escala en lugar de la verdadera estructura conjunta de los datos. Escalar las variables asegura que cada dimensión contribuya adecuadamente a la métrica de distancia, hace que el parámetro de ancho de banda sea significativo y conduce a una estimación multivariante de densidad más fiel al utilizar el KDE de sklearn.

# Escalado de variables
# ==============================================================================
scaler = StandardScaler().set_output(transform='pandas')
datos_scaled = scaler.fit_transform(datos)

Modelo KDE¶

Se utiliza validación cruzada para identificar el mejor tipo de kernel y el ancho de banda. Al disponer de pocos datos, se emplea la estrategia LeaveOneOut.

Si se dispone de una cantidad mayor de datos, es recomendable utilizar una estrategia de validación cruzada con varios folds (por ejemplo, 5 o 10).

# Validación cruzada para identificar kernel y bandwidth
# ==============================================================================
param_grid = {
    'kernel': ['gaussian', 'epanechnikov', 'exponential', 'linear'],
    'bandwidth': np.logspace(-2, 2, 20)
}

grid = GridSearchCV(
        estimator  = KernelDensity(),
        param_grid = param_grid,
        n_jobs     = -1,
        cv         = LeaveOneOut(), 
        verbose    = 0
      )
_ = grid.fit(X=datos_scaled)

resultados = (
    pd.DataFrame(grid.cv_results_)
    .loc[:, ['param_bandwidth', 'param_kernel', 'mean_test_score']]
    .sort_values(by='mean_test_score', ascending=False)
)
resultados.head(10)

# Mejores hiperparámetros por validación cruzada
# ==============================================================================
print("----------------------------------------")
print("Mejores hiperparámetros encontrados (cv)")
print("----------------------------------------")
print(grid.best_params_, ":", grid.best_score_)
modelo_kde = grid.best_estimator_

----------------------------------------
Mejores hiperparámetros encontrados (cv)
----------------------------------------
{'bandwidth': np.float64(0.18329807108324356), 'kernel': 'gaussian'} : -1.4693680944036116

Si se calcula la predicción para todo el rango de valores observado, se puede representar la función de densidad obtenida. Al tratarse de dos dimensiones, hay que crear un grid que cubra toda la superficie.

# Mapa de densidad de probabilidad
# ==============================================================================
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 4))

# Grid de valores dentro del rango observado (2 dimensiones)
x = np.linspace(min(datos['duration']) * 0.8, max(datos['duration']) * 1.1, 200)
y = np.linspace(min(datos['waiting']) * 0.8, max(datos['waiting']) * 1.1, 200)
xx, yy = np.meshgrid(x, y)
grid = pd.DataFrame(
    np.column_stack([xx.ravel(), yy.ravel()]),
    columns=['duration', 'waiting']
)
grid = scaler.transform(grid)
# Densidad de probabilidad de cada valor del grid
log_densidad_pred = modelo_kde.score_samples(grid)
densidad_pred = np.exp(log_densidad_pred).reshape(xx.shape)
ax.contour(
    xx, yy,
    densidad_pred,
    levels=12,
    alpha=0.6
)
ax.scatter(datos['duration'], datos['waiting'], color="#3182bd", alpha=0.5)
ax.set_title('Función de densidad estimada')
ax.set_xlabel('duración (duration)')
ax.set_ylabel('intervalo (waiting)');

# Representación como mapa de calor
#===============================================================================
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 4))
ax.scatter(grid['duration'], grid['waiting'], alpha=0.6, c=densidad_pred)
ax.set_title('Función de densidad estimada')
ax.set_xlabel('intervalo (waiting)')
ax.set_ylabel('duración (duration)');

# Representación 3D
#===============================================================================
fig = plt.figure(figsize=(5, 5))
ax = plt.axes(projection='3d')
ax.plot_surface(xx, yy, densidad_pred, cmap='viridis')
ax.set_xlabel('intervalo (waiting)')
ax.set_ylabel('duración (duration)')
ax.set_zlabel('densidad')
ax.set_title('KDE - Superficie 3D densidad')
plt.tight_layout();

Predicción¶

Una vez estimada la función de densidad, se puede predecir su valor para cualquier otra observación. Por ejemplo, que una erupción tenga un intervalo de 120 y una duración de 1 (esquina inferior del gráfico 3D).

new_data = np.array([120, 1]).reshape(1, 2)
log_densidad_pred = modelo_kde.score_samples(new_data)
densidad_pred = np.exp(log_densidad_pred)
densidad_pred

array([0.])

Entre los datos observados no hay ninguna erupción con estas características, de ahí que su valor de densidad sea muy bajo. De nuevo matizar que el valor obtenido es una densidad de probabilidad, no probabilidad, por lo que su rango está acotado entre [0, +$\infty$].

Este valor permite comparar cómo de verosímil es una observación frente a otra en esta misma distribución, es decir sirve para hacer comparaciones relativas.

Información de sesión¶

import session_info
session_info.show(html=False)

-----
matplotlib          3.10.8
mpl_toolkits        NA
numpy               2.2.6
pandas              2.3.3
scipy               1.15.3
seaborn             0.13.2
session_info        v1.0.1
sklearn             1.7.2
-----
IPython             9.8.0
jupyter_client      8.7.0
jupyter_core        5.9.1
-----
Python 3.13.11 | packaged by Anaconda, Inc. | (main, Dec 10 2025, 21:28:48) [GCC 14.3.0]
Linux-6.14.0-37-generic-x86_64-with-glibc2.39
-----
Session information updated at 2026-01-20 00:08

Instrucciones para citar¶

¿Cómo citar este documento?

Si utilizas este documento o alguna parte de él, te agradecemos que lo cites. ¡Muchas gracias!

Ajuste de distribuciones con kernel density estimation y Python por Joaquín Amat Rodrigo, disponible bajo una licencia Attribution-NonCommercial-ShareAlike 4.0 International (CC BY-NC-SA 4.0 DEED) en https://www.cienciadedatos.net/documentos/pystats02-ajuste-distribuciones-kde-python.html

¿Te ha gustado el artículo? Tu ayuda es importante

Tu contribución me ayudará a seguir generando contenido divulgativo gratuito. ¡Muchísimas gracias! 😊

Este documento creado por Joaquín Amat Rodrigo tiene licencia Attribution-NonCommercial-ShareAlike 4.0 International.

Se permite:

• Compartir: copiar y redistribuir el material en cualquier medio o formato.

• Adaptar: remezclar, transformar y crear a partir del material.

Bajo los siguientes términos:

• Atribución: Debes otorgar el crédito adecuado, proporcionar un enlace a la licencia e indicar si se realizaron cambios. Puedes hacerlo de cualquier manera razonable, pero no de una forma que sugiera que el licenciante te respalda o respalda tu uso.

• No-Comercial: No puedes utilizar el material para fines comerciales.

• Compartir-Igual: Si remezclas, transformas o creas a partir del material, debes distribuir tus contribuciones bajo la misma licencia que el original.