This is my first attempt for somebody coming from the frequentist camp to do Bayesian data analysis. I read a number of tutorials and few chapters from Bayesian Data Analysis by A. Gelman.
As the first more or less independent data analysis example I picked is train waiting times. I asked myself: what is the distribution of the waiting times?
The dataset was provided on on a blog and was analyzed slightly differently and outside PyMC.
My goal is to estimate the expected train wait times given those 19 data entries.
The model I built is the following:
$\mu \sim N(\hat\mu,\hat\sigma)$
$\sigma \sim |N(0,\hat\sigma)|$
$\lambda \sim \Gamma(\mu,\sigma)$
$\rho \sim Poisson(\lambda)$
where $\hat\mu$ is data mean and $\hat\sigma$ is data standard deviation multiplied by 1000.
I modeled the expected wait time as $\rho$ using the Poisson distribution. The rate parameter for this distribution is modeled using the Gamma distribution as it is conjugate distribution to Poisson distribution. The hyper-priors $\mu$ and $\sigma$ were modeled with Normal and Half-Normal distributions respectively. The standard deviation $\sigma$ was made as broad as possible to be as non-commital as possible..
I have a bunch of questions
- Is this model reasonable for the task (several possible ways to model?)?
- Did I make any beginner mistakes?
- Can the model be simplified (I tend to complicate simple things)?
- How can I verify if the posterior for the rate parameter ($\rho$) is actually fitting the data?
- How can I draw some samples from the fitted Poisson distribution to see the samples?
I can post the source code as well. In model fitting stage I do the steps for parameters $\mu$ and $\sigma$ using NUTS. Then in the second step I do Metropolis for the rate parameter $\rho$. Finally I plot the trace using the inbuilt tools.
I would be very grateful for any remarks and comments that would enable me to grasp more probabilistic programming. May be there are more classic examples that are worth experimenting with?
Here is the code I wrote in Python using PyMC3. The data file can be found here.
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()