MCMC

Exercise 1 (50 points)

Gibbs sampler example from Robert and Casella, 10.17

Suppose we have data of the number of failures (\(y_i\)) for each of 10 pumps in a nuclear plant. We also have the times (\(t_i\)) at which each pump was observed. We want to model the number of failures with a Poisson likelihood, where the expected number of failure \(\lambda_i\) differs for each pump. Since the time which we observed each pump is different, we need to scale each \(\lambda_i\) by its observed time \(t_i\). To be explicit, we assume that \(y_i\) has a Poisson distribution with rate \(\mu_i = \lambda_i t_i\).

The likelihood \(f\) is

\prod_{i=1}^{10} \text{Poisson}(\mu_i)

We let the prior \(g\) for \(\lambda\) be

\lambda \sim \text{Gamma}(\alpha_\mu, \beta_\mu)

and let the hyperprior \(h\) for \(\alpha\) to be

\alpha \sim \text{Gamma}(\alpha_\alpha, \beta_\alpha)

with \(\alpha_\alpha = 1.8\) and \(\beta_\alpha = 1.0\).

and let the hyperprior \(h\) for \(\beta\) to be

\beta \sim \text{Gamma}(\alpha_\beta, \beta_\beta)

with \(\alpha_\beta = 10.0\) and \(\beta_\beta = 1.0\).

There are 12 unknown parameters (10 \(\lambda\)s, \(\alpha\) and \(\beta\)) in this hierarchical model. Do th following using pymc3 and some plotting package.

  • Wrtie the model and run for 10,000 iterations using the No U-Turn Sampler (30 points)
  • plot the traces and distributions of the last 10% for \(\lambda_i\), \(\alpha\) and \(\beta\) (there should be 12 sets of plots) (10 points)
  • Gnnerate 1,000 samples of the number of failures \(y_i\) from the prior distribution and plot the histogram or density. That is, for each of the 10 pumps, we want to see the distribtuion of 1,000 draws of the number of failures (5 points).
  • Generate 1,000 posterior predictive samples of the number of failures \(y_i\) and plot the histogram or density. This is similar to the previous question but using draws from the posterior (5 points)

Use the following data

y = np.array([5, 1, 5, 14, 3, 19, 1, 1, 4, 22])
t = np.array([94.32, 15.72, 62.88, 125.76, 5.24, 31.44, 1.05, 1.05, 2.10, 10.48])

In [2]:

import pymc3 as pm

In [23]:

y = np.array([5, 1, 5, 14, 3, 19, 1, 1, 4, 22])
t = np.array([94.32, 15.72, 62.88, 125.76, 5.24, 31.44, 1.05, 1.05, 2.10, 10.48])

In [ ]:

Exercise 2 (50 points)

Repeat Exercise 1 using pystan.

In [ ]: