# introduccion-a-la-estadistica-bayesiana-con-pymc

Introducción a la estadística bayesiana con PyMC

Volver a la portada del curso
Notebook

Implementación en PyMC

Sección 6 · Proyecto Final

Implementación en PyMC: del dataset al reporte

Esta lección recorre, bloque por bloque, el notebook que implementa el proyecto final: una regresión lineal bayesiana sobre el dataset de diabetes de scikit-learn. Vas a comparar dos modelos, uno con un solo predictor y otro con tres, y a decidir cuál predice mejor con evidencia, no con intuición. Cada bloque es corto y se ejecuta en orden; si un paso falla, los siguientes no tienen sentido.

El modelo

Con predictores y respuesta estandarizados (media 0, desviación 1), el modelo B se escribe así:

$$y_i \sim \mathcal{N}(\mu_i, \sigma), \quad \mu_i = \alpha + \sum_{j=1}^{3} \beta_j x_{ij}$$ $$\alpha \sim \mathcal{N}(0, 1), \quad \beta_j \sim \mathcal{N}(0, 1), \quad \sigma \sim \text{HalfNormal}(1)$$

El modelo A es el caso particular con un solo predictor, bmi. Estandarizar tiene una ventaja práctica: un prior $\mathcal{N}(0, 1)$ sobre un coeficiente dice "un cambio de una desviación estándar en el predictor rara vez mueve la respuesta más de dos desviaciones estándar". Es débilmente informativo y no depende de las unidades originales.

Paso 1: carga y estandarización

load_diabetes no descarga nada: el dataset viene con scikit-learn. Tiene 442 pacientes y 10 variables (age, sex, bmi, bp, s1 a s6); la respuesta mide la progresión de la enfermedad un año después. Las variables ya vienen centradas y escaladas, pero no con desviación 1, así que se estandariza de nuevo, igual que la respuesta.

import numpy as np
import pymc as pm
import arviz as az
from sklearn.datasets import load_diabetes

datos = load_diabetes(as_frame=True)
predictores = ["bmi", "bp", "s5"]
X = datos.data[predictores]
X_std = ((X - X.mean()) / X.std()).to_numpy()
y = datos.target.to_numpy()
y_std = (y - y.mean()) / y.std()
print(X_std.shape, y_std.shape)

Paso 2: el modelo con coords y dims

Las coords dan nombre a cada dimensión; así az.summary muestra beta[bmi] en lugar de beta[0]. pm.Data guarda la matriz de diseño dentro del modelo, lo que permite cambiarla después para predecir.

coords = {"predictor": predictores, "obs": np.arange(len(y_std))}

with pm.Model(coords=coords) as modelo_b:
    X_data = pm.Data("X", X_std, dims=("obs", "predictor"))
    alpha = pm.Normal("alpha", mu=0, sigma=1)
    beta = pm.Normal("beta", mu=0, sigma=1, dims="predictor")
    sigma = pm.HalfNormal("sigma", sigma=1)
    mu = alpha + pm.math.dot(X_data, beta)
    pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y_std, dims="obs")

El modelo A se construye igual, con X_std[:, 0] como único predictor y mu = alpha + beta * x.

Paso 3: prior predictive check

Antes de ver los datos, simula respuestas desde los priors. Si generan valores de $y$ absurdos (por ejemplo, de ±50 desviaciones), los priors están mal calibrados.

with modelo_b:
    prior = pm.sample_prior_predictive(random_seed=42)

y_prior = prior.prior_predictive["y_obs"].values
print(f"Rango prior: {y_prior.min():.1f} a {y_prior.max():.1f}")

Paso 4: muestreo y log-likelihood

Fija random_seed para que el resultado sea reproducible. En Windows sin compilador de C++, usa nuts_sampler="nutpie". Para la comparación con LOO se necesita el grupo log_likelihood; con nutpie, idata_kwargs={"log_likelihood": True} no lo agrega, así que se calcula aparte con pm.compute_log_likelihood, que funciona con cualquier sampler.

with modelo_b:
    idata_b = pm.sample(
        1000, tune=1000, chains=4, cores=1,
        random_seed=42, nuts_sampler="nutpie",
    )
    pm.compute_log_likelihood(idata_b)
    pm.sample_posterior_predictive(
        idata_b, extend_inferencedata=True, random_seed=42
    )

Paso 5: diagnóstico

No interpretes nada hasta que el diagnóstico pase: r_hat ≤ 1.01, ess_bulk y ess_tail de al menos unos cientos por parámetro y cero divergencias.

resumen = az.summary(idata_b, var_names=["alpha", "beta", "sigma"])
print(resumen[["mean", "sd", "hdi_3%", "hdi_97%", "ess_bulk", "r_hat"]])

divergencias = int(idata_b.sample_stats["diverging"].sum())
print("Divergencias:", divergencias)
az.plot_trace(idata_b, var_names=["beta", "sigma"])

Paso 6: PPC y comparación con LOO

az.plot_ppc superpone la densidad de los datos observados sobre las réplicas del posterior predictivo. Después, az.compare ordena los modelos por elpd_loo: el de rango 0 es el que mejor predice datos nuevos.

az.plot_ppc(idata_b, num_pp_samples=100)

comparacion = az.compare({"A: bmi": idata_a, "B: bmi+bp+s5": idata_b}, ic="loo")
print(comparacion[["rank", "elpd_loo", "p_loo", "elpd_diff", "dse", "weight"]])
print(az.loo(idata_b))

Compara elpd_diff con dse: si la diferencia es varias veces mayor que su error estándar, la ventaja del modelo ganador es clara. Revisa también los valores de Pareto $k$ que reporta az.loo; valores por encima de 0.7 indican que la aproximación LOO no es confiable para esas observaciones.

Trampas comunes

  • Olvidar estandarizar $y$: con la respuesta en su escala original (25 a 346), un prior $\mathcal{N}(0, 1)$ para alpha y HalfNormal(1) para sigma contradice los datos.
  • Comparar sin log_likelihood: az.compare lanza un error si algún InferenceData no tiene ese grupo. Calcúlalo en ambos modelos.
  • Interpretar antes de diagnosticar: un r_hat de 1.05 o divergencias invalidan cualquier HDI que reportes.
  • Leer hdi_3% y hdi_97% como intervalo de confianza: es un intervalo de credibilidad del 94 %: contiene el parámetro con esa probabilidad posterior, dado el modelo.
  • Confundir coeficientes estandarizados con efectos en unidades: beta[bmi] se lee en desviaciones estándar de $y$ por desviación estándar de bmi.

Cierre

El orden importa tanto como el código: datos, priors verificados, muestreo, diagnóstico, PPC y, solo al final, comparación e interpretación. Si cada bloque deja una salida que puedas revisar, el reporte del proyecto se escribe casi solo a partir de esas salidas.

Recursos