Inferencia Bayesiana Aplicada: Modelos Lineales Mixtos con Bambi

Instalación de Bambi

Bambi es una capa de alto nivel: escribes una fórmula al estilo de R y obtienes un modelo de PyMC listo para muestrear. Esa comodidad depende de varias piezas que se instalan juntas. En este notebook instalas el entorno, verificas las versiones, ajustas un modelo mínimo de humo y aprendes a leer el resumen que imprime Bambi. Descarga el notebook adjunto y ejecútalo celda por celda.

Qué instala Bambi

Bambi no hace inferencia por sí mismo; coordina otras bibliotecas:

• formulae: interpreta la fórmula (y ~ x + (1|grupo)) y construye las matrices de diseño.
• PyMC: define el modelo probabilístico y ejecuta el muestreador NUTS.
• PyTensor: el motor de cálculo simbólico de PyMC; compila el modelo y calcula gradientes.
• ArviZ: guarda el resultado en un objeto InferenceData y ofrece resúmenes y gráficas (az.summary, az.plot_trace).

La documentación oficial indica que Bambi se prueba con Python 3.12 o superior. Hay dos vías de instalación:

pip install "bambi==0.17.2" "matplotlib<3.11" nutpie
# o bien, con conda
conda install -c conda-forge "bambi=0.17.2" "matplotlib<3.11" nutpie

Fija la versión. Sin ella, pip install bambi instala hoy Bambi 0.21 con PyMC 6 y ArviZ 1.x: az.summary cambia sus columnas (de hdi_3%/hdi_97% a eti89_lb/eti89_ub) y tus decimales no coincidirán con los del curso. El límite de Matplotlib evita un error de ArviZ 0.23 al importarse con Matplotlib 3.11.

Usa un entorno virtual nuevo (python -m venv .venv o un entorno de conda) para no mezclar versiones de PyMC con otros proyectos.

Paso 1: verificar versiones

Este material se probó con Bambi 0.17.2, PyMC 5.28.5, PyTensor 2.38.3, ArviZ 0.23.4 y formulae 0.6.2. La primera celda del notebook lo confirma:

import bambi as bmb
import pymc as pm
import arviz as az
import formulae
import pytensor

print("bambi", bmb.__version__)
print("pymc", pm.__version__, "| pytensor", pytensor.__version__)
print("arviz", az.__version__, "| formulae", formulae.__version__)

Si tus versiones difieren en el último dígito, los resultados serán casi idénticos. Un salto de versión mayor (por ejemplo, PyMC 6) puede cambiar nombres de argumentos.

El aviso de g++ en Windows

Al importar PyMC en Windows es común ver WARNING (pytensor.configdefaults): g++ not detected!. PyTensor intenta compilar sus operaciones en C con un compilador g++; si no lo encuentra, cae a implementaciones en Python puro y el muestreo con NUTS de PyMC se vuelve mucho más lento. No es un error: el código funciona.

Tienes dos salidas. La primera es instalar el compilador; si usas conda, el propio aviso sugiere conda install gxx. La segunda, la que usa este curso, es nutpie: un muestreador NUTS escrito en Rust que compila el modelo con Numba por defecto, así que no depende de g++. Bambi lo invoca con un solo argumento:

idata = model.fit(draws=1000, tune=1000, chains=2, cores=1,
                  random_seed=42, inference_method="nutpie")

En Bambi 0.17 el valor por defecto es inference_method="pymc" (el NUTS de PyMC). Las otras opciones documentadas son "nutpie", "blackjax", "numpyro", "vi" y "laplace".

Paso 2: modelo mínimo de humo

Una prueba de humo no busca aprender nada nuevo: solo confirma que toda la cadena (fórmula, PyMC, muestreador, ArviZ) funciona. Se simulan 100 datos con intercepto 2 y pendiente 0.5:

import numpy as np
import pandas as pd

rng = np.random.default_rng(42)
x = rng.normal(0, 1, 100)
y = 2 + 0.5 * x + rng.normal(0, 1, 100)
df = pd.DataFrame({"x": x, "y": y})

model = bmb.Model("y ~ x", df)
print(model)
idata = model.fit(draws=1000, tune=1000, chains=2, cores=1,
                  random_seed=42, inference_method="nutpie")
print(az.summary(idata).round(2))

El resumen da un intercepto con media 2.00 (HDI 94%: 1.82 a 2.19), una pendiente de 0.62 (0.38 a 0.85) y sigma de 0.99. La pendiente estimada no es exactamente 0.5 porque con 100 datos el ruido pesa, pero el valor real cae dentro del intervalo. Las columnas r_hat en 1.0 y ess_bulk por encima de 2400 indican que las cadenas convergieron.

Paso 3: leer print(model)

Antes de ajustar, imprime siempre el modelo. Para el ejemplo de humo, Bambi muestra:

       Formula: y ~ x
        Family: gaussian
          Link: mu = identity
  Observations: 100
        Priors:
    target = mu
        Common-level effects
            Intercept ~ Normal(mu: 1.9642, sigma: 2.7103)
            x ~ Normal(mu: 0.0, sigma: 3.4997)
        Auxiliary parameters
            sigma ~ HalfStudentT(nu: 4.0, sigma: 1.0818)

Cada línea responde una pregunta: Family es la distribución de la respuesta, Link la función que conecta la media con el predictor lineal y Observations cuántas filas sobrevivieron. Bajo Common-level effects están los efectos comunes con sus priors; más adelante aparecerá también Group-level effects para los efectos de grupo. Los números de los priors no son arbitrarios: Bambi genera priors débilmente informativos escalados a los datos, con una metodología similar a la del paquete rstanarm de R.

Trampas comunes

• Ignorar el aviso de g++ y esperar minutos por un modelo pequeño: cambia a inference_method="nutpie".
• Mezclar instalaciones de pip y conda en el mismo entorno: elige una sola vía.
• No fijar random_seed: sin semilla, tus decimales no coincidirán con los del notebook.
• Llamar a model.graph() sin el binario de Graphviz instalado en el sistema: el paquete de Python no basta.

Cierre

Si la prueba de humo devuelve un intercepto cercano a 2 y r_hat de 1.0, tu entorno está listo. Guarda la celda de versiones al inicio de cada notebook: es la forma más barata de diagnosticar diferencias cuando un resultado no coincide.