Nota: Mira mis anteriores artículo para una discusión práctica sobre por qué el modelado bayesiano puede ser la opción correcta para su tarea.
Este tutorial se centrará en un flujo de trabajo + guía de código para crear un modelo de regresión bayesiana en ESTANTARun lenguaje de programación probabilístico. STAN es ampliamente adoptado y se integra con el lenguaje de su elección (R, Python, shell, MATLAB, Julia, Stata). Vea el instalación guía y documentación.
Voy a utilizar el puchero para este tutorial, simplemente porque codifico en Python. Incluso si usas otro lenguaje, las prácticas bayesianas generales y la sintaxis del lenguaje STAN que analizaré aquí no varían mucho.
Para el lector más práctico, aquí hay un enlace a la computadora portátil Para este tutorial, parte de mi taller de modelado bayesiano en la Universidad Northwestern (abril de 2024).
¡Vamos a sumergirnos!
Aprendamos a construir un modelo de regresión lineal simple, el pan de cada día de cualquier estadístico, al estilo bayesiano. Suponiendo que una variable dependiente Y y covariable incógnitaPropongo el siguiente modelo simple:
Y = a + b * incógnita + ϵ
Donde ⍺ es la intersección, β es la pendiente y ϵ es un error aleatorio. Suponiendo que,
ϵ ~ Normal(0, σ)
Podemos demostrar que
Y ~Normal(a+b* INCÓGNITA, s)
Aprenderemos cómo codificar este formulario modelo en STAN.
Generar datos
Primero, generemos algunos datos falsos.
#Model Parameters
alpha = 4.0 #intercept
beta = 0.5 #slope
sigma = 1.0 #error-scale
#Generate fake data
x = 8 * np.random.rand(100)
y = alpha + beta * x
y = np.random.normal(y, scale=sigma) #noise
#visualize generated data
plt.scatter(x, y, alpha = 0.8)
Ahora que tenemos algunos datos para modelar, veamos cómo estructurarlos y pasarlos a STAN junto con las instrucciones de modelado. Esto se hace a través de modelo cadena, que normalmente contiene 4 bloques (ocasionalmente más) datos, parámetros, modeloy generado cantidadesAnalicemos cada uno de estos bloques en detalle.
Bloque de DATOS
data { //input the data to STAN
int N;
vector(N) x;
vector(N) y;
}
El datos El bloque es quizás el más simple, le dice a STAN internamente qué datos debe esperar y en qué formato. Por ejemplo, aquí pasamos:
norte: el tamaño de nuestro conjunto de datos como tipo entero. El parte declara que N≥0. (Aunque aquí es obvio que la longitud de los datos no puede ser negativa, indicar estos límites es una buena práctica estándar que puede facilitar el trabajo de STAN).
incógnita:la covariable como un vector de longitud N.
y:el dependiente como un vector de longitud N.
Ver Documentos aquí para una gama completa de tipos de datos admitidos. STAN ofrece compatibilidad con una amplia gama de tipos, como matrices, vectores, etc. Como vimos anteriormente, STAN también admite la codificación límites sobre las variables. ¡Se recomienda codificar límites! Esto conduce a modelos mejor especificados y simplifica los procesos de muestreo probabilístico que funcionan en segundo plano.
Bloque modelo
El siguiente es el modelo bloque, donde le indicamos a STAN la estructura de nuestro modelo.
//simple model block
model {
//priors
alpha ~ normal(0,10);
beta ~ normal(0,1); //model
y ~ normal(alpha + beta * x, sigma);
}
El bloque modelo también contiene un elemento importante y a menudo confuso: previo especificación. Los valores anteriores son una por excelencia parte del modelado bayesiano y debe especificarse adecuadamente para la tarea de muestreo.
Ver mi anterior artículo para una introducción al papel y la intuición detrás de los valores anteriores. Para resumir, previo es una forma funcional presupuesta para la distribución de valores de parámetros, a menudo denominada, simplemente, como creencia previa. Aunque los antecedentes no tener que exactamente fósforo La solución final, nos deben permitir muestra de ello.
En nuestro ejemplo, utilizamos valores a priori normales de media 0 con diferentes varianzas, dependiendo de cuán seguros estemos del valor medio proporcionado: 10 para alfa (muy inseguro), 1 para beta (algo seguro). Aquí, proporcioné el valor general creencia que si bien alfa puede adoptar una amplia gama de valores diferentes, la pendiente generalmente está más limitada y no tendrá una gran magnitud.
Por lo tanto, en el ejemplo anterior, la probabilidad previa para alfa es “más débil” que para beta.
A medida que los modelos se vuelven más complicados, el espacio de solución de muestreo se expande y el aporte de creencias gana importancia. De lo contrario, si no hay una intuición fuerte, es una buena práctica simplemente aportar menos creencias al modelo, es decir, utilizar una Débilmente informativo previamente y permanecer flexibles ante los datos entrantes.
La forma de y, que quizás ya hayas reconocido, es la ecuación de regresión lineal estándar.
Cantidades generadas
Por último, tenemos nuestro bloque para cantidades generadasAquí le decimos a STAN qué cantidades queremos calcular y recibir como salida.
generated quantities { //get quantities of interest from fitted model
vector(N) yhat;
vector(N) log_lik;
for (n in 1:N){
yhat(n) = normal_rng(alpha + x(n) * beta, sigma);
//generate samples from model
log_lik(n) = normal_lpdf( y(n) | alpha + x(n) * beta, sigma);
//probability of data given the model and parameters
}
}
Nota: STAN admite que los vectores se pasen directamente a ecuaciones o como iteraciones 1:N para cada elemento n. En la práctica, he descubierto que esta compatibilidad cambia con las distintas versiones de STAN, por lo que es bueno probar la declaración iterativa si la versión vectorizada no se compila.
En el ejemplo anterior:
recordar: genera muestras para y a partir de los valores de los parámetros ajustados.
registro_lik: genera probabilidad de datos dado el modelo y el valor del parámetro ajustado.
El propósito de estos valores quedará más claro cuando hablemos de la evaluación del modelo.
En total, ahora hemos especificado completamente nuestro primer modelo de regresión bayesiana simple:
model = """
data { //input the data to STAN
int N;
vector(N) x;
vector(N) y;
}
parameters {
real alpha;
real beta;
real sigma;
}model {
alpha ~ normal(0,10);
beta ~ normal(0,1);
y ~ normal(alpha + beta * x, sigma);
}generated quantities {
vector(N) yhat;
vector(N) log_lik;for (n in 1:N){
yhat(n) = normal_rng(alpha + x(n) * beta, sigma);
log_lik(n) = normal_lpdf(y(n) | alpha + x(n) * beta, sigma); }
}
"""
Todo lo que queda es compilar el modelo y ejecutar el muestreo.
#STAN takes data as a dict
data = {'N': len(x), 'x': x, 'y': y}
STAN recibe los datos de entrada en forma de diccionario. Es importante que este diccionario contenga todas las variables que le indicamos a STAN que debe esperar en el bloque de datos del modelo; de lo contrario, el modelo no se compilará.
#parameters for STAN fitting
chains = 2
samples = 1000
warmup = 10
# set seed
# Compile the model
posterior = stan.build(model, data=data, random_seed = 42)
# Train the model and generate samples
fit = posterior.sample(num_chains=chains, num_samples=samples)The .sample() method parameters control the Hamiltonian Monte Carlo (HMC) sampling process, where —
- num_cadenas: es el número de veces que repetimos el proceso de muestreo.
- num_muestras: es el número de muestras que se extraerán en cada cadena.
- calentamiento: es el número de muestras iniciales que descartamos (ya que lleva algún tiempo alcanzar la vecindad general del espacio de soluciones).
Conocer los valores correctos para estos parámetros depende tanto de la complejidad de nuestro modelo como de los recursos disponibles.
Por supuesto, los tamaños de muestreo más grandes son ideales, pero para un modelo mal especificado resultarán ser una pérdida de tiempo y de cálculo. Como anécdota, he tenido que esperar una semana para terminar de ejecutar modelos de datos grandes y descubrí que el modelo no convergía. Es importante comenzar lentamente y comprobar la coherencia del modelo antes de ejecutar un muestreo completo.
Evaluación del modelo
Las cantidades generadas se utilizan para
- evaluar la bondad del ajuste, es decir, la convergencia,
- predicciones
- Comparación de modelos
Convergencia
El primer paso para evaluar el modelo, en el marco bayesiano, es visual. Observamos los sorteos de muestra de los Montecarlo hamiltoniano (HMC) proceso de muestreo.
En términos simplistas, STAN extrae iterativamente muestras de nuestros valores de parámetros y los evalúa (HMC lo hace). forma más, pero eso está más allá de nuestro alcance actual). Para un buen ajuste, las muestras extraídas deben converger a un área general común que, idealmente, sería el mundo Óptima.
La figura de arriba muestra los muestreos para nuestro modelo en dos cadenas independientes (roja y azul).
- A la izquierda, graficamos la distribución general del valor del parámetro ajustado, es decir, posterioresEsperamos una normal distribución si el modelo y sus parámetros son bien especificado. (Por qué ¿Es eso? Bueno, una distribución normal simplemente implica que existe un cierto rango de valores de mejor ajuste para el parámetro, lo que habla en apoyo de la forma de modelo que elegimos). Además, deberíamos esperar una considerable superposición a través de cadenas SI El modelo está convergiendo hacia un óptimo.
- A la derecha, graficamos las muestras reales extraídas en cada iteración (solo para estar seguros). extra seguro). Aquí, nuevamente, deseamos ver no sólo una angosto gama pero también mucha superposición Entre los sorteos.
No todas las métricas de evaluación son visuales. Gelman et al. (1) también proponen la Barato diagnóstico que es esencialmente una medida matemática de la similitud de muestras entre cadenas. Utilizando Rhat, se puede definir un punto de corte más allá del cual las dos cadenas se consideran demasiado diferentes para converger. Sin embargo, el punto de corte es difícil de definir debido a la naturaleza iterativa del proceso y a los períodos de calentamiento variables.
La comparación visual es, por tanto, un componente crucial, independientemente de las pruebas diagnósticas.
Un pensamiento frecuentista que puede tener aquí es que, “bueno, si todo lo que tenemos son cadenas y distribuciones, ¿cuál es el valor real del parámetro?” Este es exactamente el punto. La formulación bayesiana solo trata distribucionesNO punto estimaciones con sus estadísticas de prueba difíciles de interpretar.
Dicho esto, la parte posterior todavía se puede resumir utilizando creíble intervalos como el Intervalo de alta densidad (IDH)), que incluye todos los puntos de densidad de probabilidad más alta x%.
Es importante contrastar el modelo bayesiano. creíble intervalos con frecuentista confianza intervalos.
- El intervalo creíble da una probabilidad distribución en el valores posibles Para el parámetro es decir, la probabilidad de que el parámetro asuma cada valor en algún intervalo, dados los datos.
- El intervalo de confianza se refiere a la parámetro valor como fijadoy estima en cambio la confianza de que repetido aleatorio muestreos de los datos serían fósforo.
De ahí la
El enfoque bayesiano permite que los valores de los parámetros sean fluidos y toma los datos al pie de la letra, mientras que el enfoque frecuentista exige que exista un único valor de parámetro verdadero… si tan solo tuviéramos acceso a todos los datos.
Uf. Déjalo reposar y léelo otra vez hasta que lo haga.
Otra implicación importante de utilizar intervalos creíbles, o en otras palabras, permitir que el parámetro sea variablees que las predicciones que hacemos capturan esto incertidumbre con transparenciacon un cierto porcentaje de IDH que informa la línea de mejor ajuste.
Comparación de modelos
En el marco bayesiano, la métrica de información Watanabe-Akaike (WAIC) es la opción ampliamente aceptada para la comparación de modelos. Una explicación simple de la puntuación WAIC es que estima el modelo probabilidad mientras regularizando para el número de parámetros del modelo. En palabras simples, puede explicar el sobreajuste. Este también es un atractivo importante del marco bayesiano: uno no no necesariamente necesidad Presentar un modelo validación conjunto de datos. Por lo tanto,
El modelado bayesiano ofrece una ventaja crucial cuando los datos son escasos.
La puntuación WAIC es una comparativo medida, es decir, solo tiene sentido cuando se compara con diferentes modelos que intentan explicar los mismos datos subyacentes. Por lo tanto, en la práctica, se puede seguir añadiendo más complejidad al modelo mientras el WAIC aumenta. Si en algún momento de este proceso de añadir complejidad maniática, el WAIC empieza a caer, se puede dar por terminado el asunto: cualquier otra complejidad adicional no ofrecerá una ventaja informativa para describir la distribución de datos subyacente.
Conclusión
En resumen, el bloque de modelo STAN es simplemente una cadena. Le explica a STAN lo que le vas a dar (modelo), lo que se debe encontrar (parámetros), lo que crees que está sucediendo (modelo) y lo que debería darte a cambio (cantidades generadas).
Cuando se enciende, STAN simplemente gira la manivela y emite su salida.
El verdadero desafío consiste en definir un modelo adecuado (referirse a los datos anteriores), estructurar los datos adecuadamente, preguntarle a STAN exactamente qué necesita de él y evaluar la cordura de su resultado.
Una vez que dominemos esta parte, podremos adentrarnos en el verdadero poder de STAN, donde especificar modelos cada vez más complicados se convierte en una tarea sintáctica sencilla. De hecho, en nuestro próximo tutorial haremos exactamente eso. Nos basaremos en este ejemplo de regresión simple para explorar el modelo bayesiano. Jerárquico Modelos: un estándar de la industria, de última generación, de facto… lo que sea. Veremos cómo agregar efectos aleatorios o fijos a nivel de grupo a nuestros modelos y nos maravillaremos con la facilidad de agregar complejidad mientras se mantiene la comparabilidad en el marco bayesiano.
Suscríbete si este artículo te resultó útil y para estar al tanto de más novedades.
Referencias
(1) Andrew Gelman, John B. Carlin, Hal S. Stern, David B. Dunson, Aki Vehtari y Donald B. Rubin (2013). Análisis de datos bayesianos, tercera edición.Chapman y Hall/CRC.
