pypro.mx

Introducción a PyMC

Una regresión logística clásica te da un número por coeficiente y un intervalo de confianza que casi nadie interpreta bien. La versión bayesiana te da una distribución completa para cada coeficiente y para cada predicción: puedes preguntar directamente «¿qué probabilidad hay de que la deuda aumente el riesgo de impago?» y obtener una respuesta numérica. PyMC es la biblioteca de Python que convierte esa idea en código.

Qué hace PyMC

PyMC es un lenguaje de programación probabilística. Describes el modelo como lo escribirías en papel (qué variables son desconocidas, qué distribución previa tienen y cómo generan los datos) y PyMC se encarga de calcular la distribución posterior:

$$p(\theta \mid y) \propto p(y \mid \theta)\, p(\theta)$$

donde $p(\theta)$ es el prior, $p(y \mid \theta)$ la verosimilitud y $p(\theta \mid y)$ la posterior. Para casi cualquier modelo útil esa posterior no tiene fórmula cerrada, así que PyMC la aproxima con muestras: por defecto usa NUTS, un algoritmo de Monte Carlo hamiltoniano que aprovecha gradientes calculados automáticamente por PyTensor, su motor de cómputo.

Instalación y versiones

# En la terminal:
# python -m venv .venv
# .venv\Scripts\activate          (Windows)
# source .venv/bin/activate       (macOS o Linux)
# pip install "pymc<6" "arviz<1" numba scikit-learn matplotlib

import pytensor
pytensor.config.mode = "NUMBA"   # compila con Numba, sin compilador de C

import pymc as pm
import arviz as az
print(pm.__version__, az.__version__)

El curso se probó con PyMC 5.28.5 y ArviZ 0.23.4. Fijamos arviz<1 porque ArviZ 1.x cambió los nombres de las columnas del resumen. La línea pytensor.config.mode = "NUMBA" importa: si tu equipo no tiene compilador de C (lo normal en Windows con pip), PyTensor cae a un modo Python muy lento. En nuestras pruebas, una aproximación variacional pasó de 344 s a 3.5 s con Numba, con resultados idénticos.

Anatomía de un modelo

Todo modelo vive dentro de un bloque with pm.Model(). Cada distribución que declaras ahí se registra en el modelo: las que llevan observed= forman la verosimilitud; las demás son parámetros desconocidos. Este es el modelo que usaremos en todo el curso, sobre 800 solicitudes de crédito sintéticas con tres predictores estandarizados (ingreso, razón deuda/ingreso y antigüedad laboral):

with pm.Model(coords={"pred": ["ingreso", "deuda", "antiguedad"]}) as modelo:
    X = pm.Data("X", Xs)                       # predictores estandarizados
    alpha = pm.Normal("alpha", mu=0, sigma=1.5)
    beta = pm.Normal("beta", mu=0, sigma=1, dims="pred")
    pm.Bernoulli("impago", logit_p=alpha + X @ beta, observed=y)

    idata = pm.sample(1000, tune=1000, chains=4, random_seed=42)

az.summary(idata, round_to=3)

Tres piezas merecen atención. coords y dims ponen nombre a cada coeficiente, así que el resumen dice beta[deuda] en lugar de beta[1]. pm.Data guarda los predictores en un contenedor que después puedes reemplazar para predecir clientes nuevos. Y logit_p le dice a la Bernoulli que el predictor lineal está en escala logit: es más claro que aplicar la sigmoide a mano y evita aplicarla dos veces por error (PyMC compila ambas formas a la misma log-densidad).

Qué devuelve `pm.sample`

El resultado es un objeto InferenceData de ArviZ: un contenedor con grupos (posterior, sample_stats, observed_data…) indexados por cadena y por sorteo. Con 4 cadenas de 1,000 sorteos tienes 4,000 muestras de cada parámetro. En nuestro ejemplo, az.summary da para la deuda una media de 0.625 con un HDI del 94 % entre 0.463 y 0.786, y para el ingreso −0.651 con HDI [−0.817, −0.461]. Las columnas r_hat (todas ≤ 1.002) y ess_bulk (más de 4,000) indican que las cadenas convergieron y que hay suficientes muestras efectivas; las verás a fondo en la lección de evaluación.

Como la posterior son muestras, cualquier pregunta se responde con aritmética sobre arreglos. La probabilidad de que la deuda aumente el riesgo es la fracción de muestras con coeficiente positivo:

post = idata.posterior
prob = (post["beta"].sel(pred="deuda") > 0).mean().item()
print(prob)   # 1.0: todas las muestras son positivas

Trampas comunes

No estandarizar. Con predictores en escalas muy distintas, un prior $\mathcal{N}(0, 1)$ significa cosas diferentes para cada coeficiente y el muestreador trabaja peor. Estandariza y guarda la media y la desviación para predecir después.
Olvidar observed=. Sin él, la variable es un parámetro más y PyMC simplemente muestrea del prior: el modelo «corre» pero no aprende nada de los datos.
Leer el HDI como intervalo de confianza. El HDI del 94 % es la región donde está el 94 % de la probabilidad posterior del parámetro, dado el modelo y los datos. Esa es la lectura que normalmente se le quiere dar a un intervalo de confianza, y que este no admite.
Mezclar versiones. PyMC 6 y ArviZ 1 cambian APIs y nombres de columnas; fija las versiones del curso.

Cierre

PyMC separa lo que tú decides (priors y verosimilitud) de lo que la máquina resuelve (muestrear la posterior). En el notebook siguiente construirás este modelo paso a paso y comprobarás, antes de ver los datos, si tus priors generan predicciones razonables.