Este es mi primer intento para que alguien que viene del campo frecuentista haga análisis de datos bayesianos. Leí varios tutoriales y algunos capítulos del Análisis de datos bayesianos de A. Gelman.

Como el primer ejemplo de análisis de datos más o menos independiente que elegí es el tiempo de espera del tren. 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 manera ligeramente diferente y fuera de PyMC.

Mi objetivo es estimar los tiempos de espera esperados del tren dados esos 19 datos ingresados.

El modelo que construí 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 μ es la media de datos y σ son datos de desviación estándar se multiplica por 1.000.$\hat\mu$$\hat\sigma$

Modelé el tiempo de espera esperado como utilizando la distribución de Poisson. El parámetro de velocidad para esta distribución se modela utilizando la distribución Gamma, ya que es la distribución conjugada a la distribución de Poisson. Los hiper-previos μ y σ fueron modelados con distribuciones Normal y Media Normal, respectivamente. La desviación estándar σ se hizo lo más amplia posible para que fuera lo menos comprometida posible.$\rho$$\mu$$\sigma$$\sigma$

Tengo muchas preguntas

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

Los posteriores después de 5000 pasos de Metropolis se ven así:

Puedo publicar el código fuente también. En la etapa de ajuste del modelo, realizo los pasos para los parámetros y σ utilizando NUTS. Luego, en el segundo paso, hago Metrópolis para el parámetro de velocidad ρ . Finalmente trazo el rastro usando las herramientas incorporadas.$\mu$$\sigma$$\rho$

Estaría muy agradecido por cualquier comentario y comentario que me permitiera captar una programación más probabilística. ¿Puede haber ejemplos más clásicos con los que valga 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()

Primero te diré lo que haría y luego responderé las preguntas específicas que tengas.

Lo que haría (al menos inicialmente)

Esto es lo que deduzco de su publicación: tiene tiempos de espera de capacitación para 19 observaciones y está interesado en inferencia sobre el tiempo de espera esperado.

$W_i$$i=1,\ldots,19$$i$$W_i\in\mathbb{R}^+$

Hay varios supuestos posibles del modelo que podrían usarse y con 19 observaciones puede ser difícil determinar qué modelo es más razonable. Algunos ejemplos son log-normal, gamma, exponencial, Weibull.

$Y_i=\log(W_i)$

$E[W_i] = e^{\mu+\sigma^/2}$$\mu$$\sigma^2$$e^{\mu+\sigma^/2}$

Contestando tus preguntas

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

$\lambda$$\lambda$

• ¿Cometí algún error de principiante?

Ver comentario anterior.

$\lambda$

Su prior no debe depender de los datos.

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

Si y deberia. Vea mi enfoque de modelado.

• ¿Cómo puedo verificar si la parte posterior del parámetro de velocidad ( ρ$\rho$

$\rho$$\lambda$

• ¿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 su MCMC, conecte los valores de los parámetros para esa iteración y tome una muestra.