Inferencia estadística avanzada en Python: bootstrap, máxima verosimilitud y modelos lineales generalizados

Cuando los supuestos clásicos (normalidad, varianza conocida, muestras grandes) no se cumplen, la estadística moderna recurre a dos herramientas poderosas: la estimación por máxima verosimilitud (MLE) y los métodos de remuestreo como el bootstrap. En este artículo las implementamos en Python con SciPy y statsmodels, y cerramos con un modelo lineal generalizado (GLM) para datos de conteo. Es un recorrido de nivel intermedio-avanzado, orientado a quien ya domina numpy y pandas.

1. Máxima verosimilitud (MLE)

Dada una muestra \(x_1, \dots, x_n\) i.i.d. de una densidad \(f(x\mid\theta)\), la función de verosimilitud es \(L(\theta) = \prod_{i=1}^{n} f(x_i\mid\theta)\). Por estabilidad numérica se maximiza su logaritmo:

\[\hat{\theta}_{\text{MLE}} = \arg\max_{\theta} \; \ell(\theta), \qquad \ell(\theta) = \sum_{i=1}^{n}\log f(x_i\mid\theta)\]

Bajo condiciones de regularidad, el estimador MLE es consistente y asintóticamente normal, con varianza dada por la inversa de la información de Fisher \(\mathcal{I}(\theta)\):

\[\sqrt{n}\,(\hat{\theta} - \theta) \xrightarrow{d} \mathcal{N}\!\left(0,\; \mathcal{I}(\theta)^{-1}\right)\]

Ajustemos los parámetros de una distribución Gamma a datos simulados, maximizando la log-verosimilitud con scipy.optimize:

import numpy as np
from scipy import stats, optimize

rng = np.random.default_rng(42)
data = rng.gamma(shape=2.0, scale=1.5, size=500)  # verdad: k=2.0, θ=1.5

def neg_log_likelihood(params, x):
    shape, scale = params
    if shape <= 0 or scale <= 0:
        return np.inf
    return -np.sum(stats.gamma.logpdf(x, a=shape, scale=scale))

result = optimize.minimize(
    neg_log_likelihood,
    x0=[1.0, 1.0],
    args=(data,),
    method="Nelder-Mead",
)
shape_hat, scale_hat = result.x
print(f"shape (k) estimado = {shape_hat:.3f}")
print(f"scale (θ) estimado = {scale_hat:.3f}")

Salida típica:

shape (k) estimado = 2.034
scale (θ) estimado = 1.472

Los valores recuperan la verdad poblacional (\(k=2.0\), \(\theta=1.5\)) dentro del error muestral.

2. Bootstrap: intervalos de confianza sin fórmulas

El bootstrap de Efron estima la distribución muestral de un estadístico remuestreando con reemplazo de los propios datos. Es ideal cuando no existe una fórmula cerrada para el error estándar (por ejemplo, para la mediana o para un coeficiente de correlación).

El procedimiento para \(B\) réplicas es:

1. Muestrear \(n\) observaciones con reemplazo de la muestra original.
2. Calcular el estadístico \(\hat{\theta}^{*}_b\) sobre cada réplica.
3. Usar la distribución de \(\{\hat{\theta}^{*}_b\}_{b=1}^{B}\) para construir el intervalo (método percentil).

import numpy as np

rng = np.random.default_rng(7)
sample = rng.lognormal(mean=1.0, sigma=0.5, size=200)

def bootstrap_ci(x, statistic, n_boot=10_000, alpha=0.05, seed=0):
    rng = np.random.default_rng(seed)
    n = len(x)
    estimates = np.empty(n_boot)
    for b in range(n_boot):
        resample = rng.choice(x, size=n, replace=True)
        estimates[b] = statistic(resample)
    lower = np.percentile(estimates, 100 * alpha / 2)
    upper = np.percentile(estimates, 100 * (1 - alpha / 2))
    return statistic(x), (lower, upper), estimates

theta_hat, (lo, hi), dist = bootstrap_ci(sample, np.median)
print(f"Mediana = {theta_hat:.3f}")
print(f"IC 95% (percentil) = [{lo:.3f}, {hi:.3f}]")

Mediana = 2.701
IC 95% (percentil) = [2.447, 2.985]

La distribución bootstrap del estadístico (abajo) suele ser aproximadamente normal por el Teorema del Límite Central, incluso partiendo de datos asimétricos (aquí, log-normales):

3. Modelos lineales generalizados (GLM)

Los GLM extienden la regresión lineal a respuestas no normales (conteos, binarias, proporciones) mediante una función de enlace \(g(\cdot)\) que conecta la media \(\mu\) con el predictor lineal:

\[g(\mu) = X\beta, \qquad \mathbb{E}[y\mid X] = \mu = g^{-1}(X\beta)\]

Para datos de conteo se usa la regresión de Poisson, con enlace logarítmico \(g(\mu) = \log(\mu)\), de modo que \(\mu = e^{X\beta}\) garantiza una media positiva. Ajustémoslo con statsmodels:

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

rng = np.random.default_rng(123)
n = 300
horas_estudio = rng.uniform(0, 10, n)
# Verdad: log(μ) = 0.3 + 0.25 * horas
mu = np.exp(0.3 + 0.25 * horas_estudio)
ejercicios = rng.poisson(mu)

df = pd.DataFrame({"ejercicios": ejercicios, "horas": horas_estudio})

modelo = smf.glm(
    formula="ejercicios ~ horas",
    data=df,
    family=sm.families.Poisson(),
).fit()

print(modelo.summary())

La salida de statsmodels incluye coeficientes, errores estándar, estadísticos de Wald e intervalos de confianza:

                 Generalized Linear Model Regression Results
==============================================================================
Dep. Variable:             ejercicios   No. Observations:                  300
Model:                            GLM   Df Residuals:                      298
Model Family:                 Poisson   Df Model:                            1
Link Function:                    Log   Scale:                          1.0000
Method:                          IRLS   Log-Likelihood:                -612.34
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept      0.3118      0.055      5.669      0.000       0.204       0.420
horas          0.2487      0.008     30.942      0.000       0.233       0.264
==============================================================================

El coeficiente de horas (\(\hat{\beta}_1 \approx 0.249\)) se interpreta de forma multiplicativa: \(e^{0.249} \approx 1.28\), es decir, cada hora adicional de estudio multiplica por ~1.28 el número esperado de ejercicios resueltos. Ambos estimadores recuperan la verdad poblacional (0.3 y 0.25) con \(p < 0.001\).

Diagnóstico y sobredispersión

Un supuesto clave de Poisson es que media y varianza son iguales (\(\mathbb{E}[y]=\text{Var}[y]=\mu\)). Cuando la varianza es mayor (sobredispersión), conviene comparar con un modelo binomial negativo y revisar el estadístico de Pearson:

pearson_chi2 = modelo.pearson_chi2
dispersion = pearson_chi2 / modelo.df_resid
print(f"Dispersión (≈1 ideal) = {dispersion:.3f}")

Un valor cercano a 1 confirma que Poisson es adecuado; valores claramente superiores a 1 sugieren cambiar de familia o añadir un término de dispersión.

4. Conclusión

MLE, bootstrap y GLM forman un núcleo práctico de la estadística computacional: el primero da estimadores óptimos bajo un modelo paramétrico, el segundo cuantifica incertidumbre sin fórmulas analíticas, y el tercero generaliza la regresión a respuestas de cualquier familia exponencial. La combinación de SciPy (optimización y distribuciones) y statsmodels (modelos con inferencia completa) cubre la mayoría de los flujos de análisis sin salir de Python.