Guía de Trabajo Práctico 06

Guía de Trabajo Práctico 06#

Materia: Aprendizaje Profundo Basado en la Física (optativa del Instituto Balseiro)
Docente: José I. Robledo
Edición: abril-mayo 2026

Inferencia Bayesiana basada en simulaciones#

En esta guía vas a trabajar con problemas de simulation-based inference (SBI), donde el modelo físico o ingenieril está disponible como simulador pero la verosimilitud puede ser difícil, costosa o directamente inaccesible.

La secuencia queda organizada así:

un problema simple para comparar inferencia bayesiana tradicional con ABC;
un simulador de conteos radiactivos para implementar ABC por rechazo;
el mismo detector resuelto con SNPE, estadísticas resumen, embeddings y un posterior predictive check;
un oscilador de Duffing como ejemplo de simulador físico no lineal y construcción de estadísticas resumen;
una comparación entre SNPE, SNLE y SNRE sobre el problema resumido de Duffing.

Usaremos la expresión estadísticas resumen como traducción de summary statistics. En algunos bloques aparece también el término en inglés porque es habitual en la documentación de sbi.

La idea es ir desde un caso controlado, donde conocemos la posterior exacta, hacia problemas donde necesitamos compresión de datos, aproximaciones neuronales y diagnósticos predictivos.

Si el entrenamiento de sbi tarda demasiado, reducí primero num_simulations y max_num_epochs, y recién después aumentalos para mejorar la calidad del posterior.

import numpy as np
import matplotlib.pyplot as plt
import torch
import torch.nn as nn
from scipy import stats
from scipy.integrate import solve_ivp
from scipy.signal import find_peaks

from sbi.inference import SNPE, infer
from sbi.analysis import pairplot
from sbi.neural_nets import posterior_nn
from sbi.utils import BoxUniform

seed = 42
rng = np.random.default_rng(seed)
torch.manual_seed(seed)
np.random.seed(seed)

plt.rcParams["figure.figsize"] = (7, 4)
plt.rcParams["axes.grid"] = True

device = "cuda" if torch.cuda.is_available() else "cpu"
print("device:", device)

Ejercicio 1 — Moneda sesgada: Bayes exacto vs ABC#

En este primer problema vamos a inferir la probabilidad $p$ de obtener cara al lanzar una moneda. Es un ejemplo ideal para comparar dos enfoques:

inferencia bayesiana tradicional, donde podemos escribir la posterior exacta porque el modelo Bernoulli con prior Beta es conjugado;
ABC por rechazo, donde evitamos usar la verosimilitud explícita y aceptamos parámetros cuyas simulaciones se parecen a los datos observados.

Trabajá con n_trials = 50 lanzamientos y una prior uniforme, equivalente a $\mathrm{Beta}(1,1)$.

Generar una secuencia de observaciones con un valor verdadero p_true elegido por ustedes.
Calcular la posterior exacta $\mathrm{Beta}(\alpha,\beta)$ luego de observar los datos.
Graficar la posterior exacta para distintos tamaños de muestra: 0, 1, 2, 5, 10, 20, 50 lanzamientos.
Implementar un ABC por rechazo usando como simulador una moneda con parámetro p y como distancia la discrepancia entre el número de caras simulado y observado.
Repetir ABC con al menos dos tolerancias distintas o aceptando distintos porcentajes de muestras.
Comparar el histograma de muestras aceptadas por ABC con la posterior exacta.
Discutir por qué, en este problema, el número total de caras funciona como estadística suficiente.

p_true = 0.63
n_trials = 50
coin_data = stats.bernoulli.rvs(p_true, size=n_trials, random_state=seed)
ps = np.linspace(0.001, 0.999, 500)
heads_obs = int(coin_data.sum())

# 1. Calcular alpha_post y beta_post para distintos prefijos de coin_data.
# 2. Graficar las distribuciones Beta exactas.
# 3. Implementar un simulador ABC para el numero de caras.
# 4. Muestrear muchos valores de p desde la prior y aceptar los que reproduzcan mejor heads_obs.
# 5. Comparar las muestras aceptadas contra la posterior exacta.

# Sugerencia:
# alpha_post = 1 + heads_obs
# beta_post = 1 + n_trials - heads_obs
# exact_pdf = stats.beta.pdf(ps, alpha_post, beta_post)

Ejercicio 2 — ABC por rechazo en conteos radiactivos#

Ahora consideremos un detector que registra conteos de una fuente radiactiva en distintas ventanas temporales. El número esperado de conteos en el tiempo $t_i$ se modela como

\[ \mu_i(\lambda, \epsilon, b) = N_0\,\epsilon\,\exp(-\lambda t_i) + b, \]

donde $\lambda$ es la tasa de decaimiento, $\epsilon$ es la eficiencia del detector y $b$ es un fondo aproximadamente constante. La observación se genera como

\[ x_i \sim \mathrm{Poisson}(\mu_i). \]

Definir una prior uniforme sobre $\theta=(\lambda,\epsilon,b)$ usando rangos físicamente razonables.
Implementar un simulador vectorizado que acepte un tensor theta de forma (n, 3) y devuelva conteos de forma (n, n_times).
Generar una observación sintética x_obs a partir de un theta_true=(0.42, 0.78, 6.0).
Implementar ABC por rechazo:
- muestrear muchos parámetros del prior;
- simular datos para cada parámetro;
- calcular una distancia entre simulación y observación;
- aceptar el 1%, 2% o 5% más cercano.
Repetir el punto anterior usando estadísticas resumen en lugar de las observaciones completas. Por ejemplo: conteo total, razón temprano/tardío y pendiente de una regresión lineal de log(counts + 1) contra el tiempo.
Graficar las muestras aceptadas en pares de parámetros e indicar dónde está theta_true.

times_counts = torch.linspace(0.5, 6.0, 12)
N0 = 350.0
theta_true_counts = torch.tensor([[0.42, 0.78, 6.0]])

prior_counts = BoxUniform(
    low=torch.tensor([0.05, 0.35, 0.0]),
    high=torch.tensor([1.20, 0.98, 18.0]),
)


def radioactive_simulator(theta):
    """Simulador vectorizado de conteos Poisson.

    Parametros
    ----------
    theta : torch.Tensor, shape (n, 3)
        Columnas: lambda, eficiencia, fondo.
    """
    # Completar: construir mean_counts con broadcasting y devolver torch.poisson(mean_counts).
    # Recordar que la salida debe tener shape (n, len(times_counts)).
    raise NotImplementedError


def count_summaries(x):
    """Estadisticas resumen para una matriz de conteos."""
    # Completar con 3 a 5 estadisticas informativas y escalarlas si hace falta.
    # Sugerencia: total, media temprana, media tardia, razon temprana/tardia, pendiente log-lineal.
    raise NotImplementedError

Ejercicio 3 — SNPE para el detector: estadísticas resumen, embeddings y PPC#

Ahora reutilizá el simulador del detector radiactivo para entrenar una posterior neuronal con SNPE.

Generar un conjunto de entrenamiento con pares (theta, x) simulados desde prior_counts.
Entrenar SNPE usando las estadísticas resumen de count_summaries(x).
Definir una embedding_net pequeña y entrenar otro SNPE directamente sobre los conteos completos.
Comparar visualmente los posteriores obtenidos con estadísticas resumen manuales y con embedding aprendido.
Calcular media posterior, desvío estándar e intervalo creíble central del 90% para cada parámetro.
Realizar un posterior predictive check (PPC) con al menos 100 muestras posteriores:
- muestrear parámetros desde la posterior;
- volver a simular conteos;
- graficar la observación, la mediana predictiva y una banda central del 90%.
Discutir qué gana y qué pierde cada enfoque: interpretabilidad de las estadísticas resumen manuales frente a flexibilidad del embedding aprendido.

Ejercicio 4 — Oscilador de Duffing#

En muchos problemas de física e ingeniería, una simulación devuelve una trayectoria temporal larga y no un vector corto. Para practicar ese caso vamos a usar un oscilador de Duffing amortiguado y forzado.

El modelo para el desplazamiento $x(t)$ es

\[ \ddot{x} + 2\gamma \dot{x} + \omega_0^2 x + \beta x^3 = F\cos(\omega_d t), \]

donde vamos a inferir $$ \theta = (\gamma, \beta, F, \sigma), $$ con $\sigma$ representando ruido gaussiano de medición.

El objetivo de este ejercicio no es todavía entrenar un método de sbi, sino pensar cómo resumir una trayectoria de manera razonable.

Ejecutar el simulador provisto para varios valores de theta y visualizar las trayectorias $x(t)$.
Generar una observación sintética con un theta_true_duffing conocido.
Construir una función summarize_duffing(x) con un conjunto corto de estadísticas resumen informativas. Por ejemplo:
- media y desvío estándar;
- amplitud máxima;
- RMS de la parte estacionaria;
- frecuencia dominante estimada con la FFT;
- cantidad de máximos detectados con find_peaks.
Comparar la dimensión de la trayectoria completa con la dimensión del resumen.
Explicar brevemente qué parámetros deberían quedar mejor informados por cada estadística y dónde esperan degeneraciones.

t_duffing = np.linspace(0.0, 50.0, 500)
x0_duffing = np.array([0.15, 0.0])
omega0_duffing = 1.0
omega_drive_duffing = 0.85


def solve_duffing(theta, t_eval=t_duffing, x0=x0_duffing):
    """Simular posicion observada de un oscilador de Duffing.

    Parametros
    ----------
    theta : array-like, shape (4,)
        gamma, beta, F, sigma.
    """
    gamma, beta, force, sigma = [float(v) for v in theta]

    def rhs(t, y):
        x, v = y
        drive = force * np.cos(omega_drive_duffing * t)
        dxdt = v
        dvdt = drive - 2.0 * gamma * v - omega0_duffing**2 * x - beta * x**3
        return [dxdt, dvdt]

    sol = solve_ivp(
        rhs,
        (t_eval[0], t_eval[-1]),
        x0,
        t_eval=t_eval,
        rtol=1e-6,
        atol=1e-8,
        max_step=0.1,
    )
    clean_position = sol.y[0].astype(np.float32)
    noisy_position = clean_position + rng.normal(0.0, sigma, size=clean_position.shape).astype(np.float32)
    return noisy_position


def plot_duffing_trace(position, title="Oscilador de Duffing"):
    plt.plot(t_duffing, position, label="posicion observada")
    plt.plot(
        t_duffing,
        0.15 * np.cos(omega_drive_duffing * t_duffing),
        "--",
        alpha=0.7,
        label="referencia de fase",
    )
    plt.xlabel("tiempo")
    plt.ylabel("x(t)")
    plt.title(title)
    plt.legend()
    plt.tight_layout()
    plt.show()


def summarize_duffing(position):
    """Convertir una trayectoria x(t) en un vector corto de estadisticas resumen."""
    # Completar con 5 o 6 estadisticas simples e informativas.
    # Sugerencia: media, std, amplitud maxima, RMS estacionario,
    # frecuencia dominante y numero de picos.
    raise NotImplementedError


def duffing_summary_simulator(theta):
    """Wrapper compatible con sbi: theta tiene shape (n, 4)."""
    summaries = []
    for row in theta.detach().cpu().numpy():
        position = solve_duffing(row)
        summaries.append(summarize_duffing(position))
    return torch.as_tensor(np.stack(summaries), dtype=torch.float32)

Ejercicio 5 — Comparar SNPE, SNLE y SNRE#

Con el simulador de Duffing resumido, vamos a comparar tres familias de métodos disponibles en sbi:

SNPE: aproxima directamente $p(\theta \mid x)$;
SNLE o NLE: aproxima $p(x \mid \theta)$ y luego construye el posterior;
SNRE o NRE: aproxima razones de verosimilitud mediante clasificación.

Definir una prior para $\gamma$, $\beta$, $F$ y $\sigma$.
Entrenar SNPE con SNPE(...) o con infer(..., method="SNPE").
Entrenar SNLE con infer(..., method="SNLE") usando el mismo número de simulaciones.
Entrenar SNRE con infer(..., method="SNRE") usando el mismo número de simulaciones.
Muestrear los tres posteriores condicionados a x_obs_duffing_summary.
Comparar los tres resultados con pairplot o con gráficos de dispersión por pares.
Registrar el tiempo de entrenamiento y el tiempo de muestreo posterior. ¿Qué método fue más barato? ¿Cuál produjo muestras más estables?
Repetir con menos simulaciones y describir cuál de los métodos se degrada primero.
Interpretar físicamente las correlaciones posteriores: ¿aparecen compensaciones entre amortiguamiento y fuerza externa? ¿La rigidez no lineal queda bien identificada?

prior_duffing = BoxUniform(
    low=torch.tensor([0.02, 0.00, 0.05, 0.005]),
    high=torch.tensor([0.22, 4.00, 0.65, 0.080]),
)