Modelización bayesiana de los tiempos de espera de los trenes: La definición del modelo

Question

Este es mi primer intento para alguien que viene del campo frecuentista de hacer análisis de datos bayesianos. He leído varios tutoriales y algunos capítulos de Bayesian Data Analysis de A. Gelman.

El primer ejemplo de análisis de datos más o menos independiente que he elegido es el de los tiempos de espera de los trenes. Me pregunté: ¿cuál es la distribución de los tiempos de espera?

El conjunto de datos se proporcionó en un blog y se analizó de forma ligeramente diferente y fuera de PyMC.

Mi objetivo es estimar los tiempos de espera esperados del tren a partir de esos 19 datos.

El modelo que he construido es el siguiente:

$\mu \sim N(\hat\mu,\hat\sigma)$

$\sigma \sim |N(0,\hat\sigma)|$

$\lambda \sim \Gamma(\mu,\sigma)$

$\rho \sim Poisson(\lambda)$

donde $\hat\mu$ es la media de los datos y $\hat\sigma$ es la desviación estándar de los datos multiplicada por 1000.

He modelado el tiempo de espera previsto como $\rho$ utilizando la distribución de Poisson. El parámetro de tasa para esta distribución se modela utilizando la distribución Gamma, ya que es una distribución conjugada a la distribución Poisson. Los hiperprecios $\mu$ y $\sigma$ se modelaron con las distribuciones Normal y Seminormal, respectivamente. La desviación estándar $\sigma$ se hizo lo más amplio posible para ser lo menos comprometido posible..

Tengo un montón de preguntas

• ¿Es este modelo razonable para la tarea (varias formas posibles de modelar?)?
• ¿He cometido algún error de principiante?
• ¿Se puede simplificar el modelo (tiendo a complicar las cosas sencillas)?
• ¿Cómo puedo verificar si la posterior para el parámetro de la tasa ( $\rho$ ) se ajusta realmente a los datos?
• ¿Cómo puedo extraer algunas muestras de la distribución de Poisson ajustada para ver las muestras?

Los posteriors después de 5000 pasos de Metrópolis se ven así:

También puedo publicar el código fuente. En la etapa de ajuste del modelo hago los pasos para los parámetros $\mu$ y $\sigma$ usando NUTS. Luego, en el segundo paso, hago Metrópolis para el parámetro de la tasa $\rho$ . Por último, trazo la huella utilizando las herramientas incorporadas.

Estaría muy agradecido por cualquier observación y comentario que me permita comprender mejor la programación probabilística. ¿Quizás haya más ejemplos clásicos con los que merezca la pena experimentar?

---

Aquí está el código que escribí en Python usando PyMC3. El archivo de datos se puede encontrar aquí .

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import pymc3

from scipy import optimize

from pylab import figure, axes, title, show

from pymc3.distributions import Normal, HalfNormal, Poisson, Gamma, Exponential
from pymc3 import find_MAP
from pymc3 import Metropolis, NUTS, sample
from pymc3 import summary, traceplot

df = pd.read_csv( 'train_wait.csv' )

diff_mean = np.mean( df["diff"] )
diff_std = 1000*np.std( df["diff"] )

model = pymc3.Model()

with model:
    # unknown model parameters
    mu = Normal('mu',mu=diff_mean,sd=diff_std)
    sd = HalfNormal('sd',sd=diff_std)

    # unknown model parameter of interest
    rate = Gamma( 'rate', mu=mu, sd=sd )

    # observed
    diff = Poisson( 'diff', rate, observed=df["diff"] )

with model:
    step1 = NUTS([mu,sd])
    step2 = Metropolis([rate])
    trace = sample( 5000, step=[step1,step2] )

plt.figure()
traceplot(trace)
plt.savefig("rate.pdf")
plt.show()
plt.close()

Answer 1

Primero te diré lo que yo haría y luego responderé a las preguntas concretas que tenías.

Lo que yo haría (al menos inicialmente)

Esto es lo que deduzco de tu post, tienes tiempos de espera de entrenamiento para 19 observaciones y estás interesado en la inferencia sobre el tiempo de espera esperado.

Definiré $W_i$ para $i=1,\ldots,19$ el tiempo de espera del tren $i$ . No veo ninguna razón para que estos tiempos de espera sean enteros, así que supondré que son cantidades continuas positivas, es decir $W_i\in\mathbb{R}^+$ . Supongo que todos los tiempos de espera se cumplen realmente.

Hay varias hipótesis de modelo posibles que podrían utilizarse y con 19 observaciones puede ser difícil determinar qué modelo es más razonable. Algunos ejemplos son log-normal, gamma, exponencial, Weibull.

Como primer modelo, sugeriría modelar $Y_i=\log(W_i)$ y luego asumir $$ Y_i \stackrel{ind}{\sim} N(\mu,\sigma^2). $$ Con esta elección, se tiene acceso a la riqueza de la teoría normal que existe, por ejemplo, una prioridad conjugada. La prioridad conjugada es una distribución normal-inversa-gamma, es decir $$ \mu|\sigma^2 \sim N(m,\sigma^2 C) \quad \sigma^2 \sim IG(a,b) $$ donde $IG$ es la distribución gamma inversa. Como alternativa, se puede utilizar la prioridad por defecto $p(\mu,\sigma^2)\propto 1/\sigma^2$ en cuyo caso la posterior es también una distribución normal-inversa-gamma.

Desde $E[W_i] = e^{\mu+\sigma^/2}$ podemos responder a preguntas sobre el tiempo de espera esperado extrayendo muestras conjuntas de $\mu$ y $\sigma^2$ a partir de su distribución posterior, que es una distribución normal-inversa-gamma, y luego calcular $e^{\mu+\sigma^/2}$ para cada una de estas muestras. Esto se muestrea a partir de la posterior para el tiempo de espera esperado.

Responder a sus preguntas

• ¿Es este modelo razonable para la tarea (varias formas posibles de modelar?)?

Un Poisson no parece apropiado para datos que podrían tener valores no enteros. Sólo tiene un único $\lambda$ y, por lo tanto, no puede aprender los parámetros de la distribución gamma que ha asignado a $\lambda$ . Otra forma de decir esto es que se ha construido un modelo jerárquico, pero no hay estructura jerárquica en los datos.

• ¿He cometido algún error de principiante?

Véase el comentario anterior.

Además, será de gran ayuda si tus matemáticas y tu código coinciden, por ejemplo, dónde está $\lambda$ en sus resultados de MCMC? ¿qué es sd y rate en su código?

Su prioridad no debe depender de los datos.

• ¿Se puede simplificar el modelo (tiendo a complicar las cosas sencillas)?

Sí y debería. Vea mi enfoque de modelado.

• ¿Cómo puedo verificar si el posterior para el parámetro de la tasa ( $\rho$ ) se ajusta realmente a los datos?

¿No es $\rho$ ¿se supone que son sus datos? ¿Quiere decir que $\lambda$ ? Una cosa que hay que comprobar es que el tiempo medio de espera de la muestra tiene sentido en relación con su distribución posterior sobre el tiempo medio de espera. A no ser que tengas un prior extraño, la media muestral debería estar cerca del pico de la distribución posterior.

• ¿Cómo puedo extraer algunas muestras de la distribución de Poisson ajustada para ver las muestras?

Creo que quieres una distribución predictiva posterior. Para cada iteración en tu MCMC, introduces los valores de los parámetros para esa iteración y tomas una muestra.