Change-point analysis with MCMC
Data: Time series of mining accidents
Techniques: Bayesian analysis, MCMC
Mining disasters
This project uses data from a time series of recorded coal mining disasters in the UK from 1851 to 1962. Sometime during this period, new legislation was passed that changed safety regulations. The number of disasters is thought to have gone down because of these changes, but we don’t know for sure when this change-point (switchpoint) occurred.
Here I go through a tutorial to estimate when the change-point occurred, and the rate of disasters before and after the change-point, using PYMC3. There will be a lot of annotations to explain what the code is doing.
%matplotlib inline
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from pymc3 import DiscreteUniform, Poisson, Exponential
from pymc3.math import switch
from pymc3 import Metropolis, NUTS, sample, Model, traceplot
from pymc3 import summary
sns.set_style("whitegrid")
sns.set_context("poster")
# The data
disaster_data = np.ma.masked_values([4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
2, 2, 3, 4, 2, 1, 3, -999, 2, 1, 1, 1, 1, 3, 0, 0,
1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
3, 3, 1, -999, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1], value=-999)
year = np.arange(1851, 1962)
plt.plot(year, disaster_data, 'o', markersize=8);
plt.ylabel("Disaster count")
plt.xlabel("Year")
The plot of disasters over time seems to show that the change-point occurred sometime in the early 1900s.
Specifying the model
We can model the occurrences of disasters with a Poisson, with an early rate for the early part of the time series, and a later (smaller) rate for the later part.
For the year in which the switchpoint occurred, we are just using a uniform distribution from the beginning to the end of the period (this is a rather “diffuse” prior that does not give the model much information on when the switchpoint is thought to have occurred).
# Specifying the model
with Model() as disaster_model:
# Discrete uniform prior for the switchpoint
switchpoint = DiscreteUniform('switchpoint', lower=year.min(), upper=year.max(), testval=1900)
# Exponential priors for pre- and post-switch rates
early_rate = Exponential('early_rate', 1)
late_rate = Exponential('late_rate', 1)
# <switch> in PYMC3 works like an if statement- it uses the first argument (switchpoint) to switch between the next two arguments.
# Allocate appropriate Poisson rates to years before and after the current
rate = switch(switchpoint >= year, early_rate, late_rate)
disasters = Poisson('disasters', rate, observed=disaster_data)
Applied log-transform to early_rate and added transformed early_rate_log_ to model.
Applied log-transform to late_rate and added transformed late_rate_log_ to model.
The Model() object is a container for the model random variables.
All statements indented after the ‘with’ statement specify the components of the model.
We have several stochastic variables here, with Poisson, discrete-uniform, and exponential priors. They are stochastic because their values are partly random.
We also have a deterministic random variable (rate). “Deterministic” (as opposed to stochastic) means that its value is completely determined by its parents’ values.
Obtain posterior estimates
with disaster_model:
trace = sample(10000) # draw 10000 posterior samples.
Assigned Metropolis to switchpoint
Assigned NUTS to early_rate_log_
Assigned NUTS to late_rate_log_
Assigned Metropolis to disasters_missing
[-----------------100%-----------------] 10000 of 10000 complete in 6.1 sec
The ‘sample’ function runs the step method (a particular MCMC sampling algorithm). These can be assigned manually or assigned automatically by PYMC3.
-Binary variables will be assigned to BinaryMetropolis
-Discrete variables will be assigned to Metropolis
-Continuous variables will be assigned to NUTS
It returns a trace object, containing the samples collected.
Posterior analysis
Let’s look at the output from the sampling.
The left column: a smoothed histogram (using kernel density estimation) of the marginal posteriors of each stochastic random variable. The right column: the samples of the Markov chain plotted in order.
# Look at the posterior plot
traceplot(trace);
There is about a 10 year span that’s credible for our switchpoint, though it looks like most of the probability mass is over a 5 year span around the early 1890s- this is our interval estimate of when the switchpoint occurred.
Posterior statistics
We can also get an output of some common posterior statistics.
summary(trace)
switchpoint:
Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
1889.726 2.428 0.097 [1884.000, 1894.000]
Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
1886.000 1888.000 1890.000 1891.000 1896.000
disasters_missing:
Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
2.075 1.819 0.081 [0.000, 6.000]
0.921 0.939 0.026 [0.000, 3.000]
Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
0.000 1.000 2.000 3.000 6.000
0.000 0.000 1.000 1.000 3.000
early_rate:
Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
3.093 0.287 0.005 [2.550, 3.651]
Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
2.561 2.895 3.085 3.280 3.667
late_rate:
Mean SD MC Error 95% HPD interval
-------------------------------------------------------------------
0.930 0.119 0.002 [0.696, 1.158]
Posterior quantiles:
2.5 25 50 75 97.5
|--------------|==============|==============|--------------|
0.715 0.848 0.924 1.005 1.180
We estimate that the switchpoint occurred between 1884 and 1894- this is the highest posterior density interval.
The early rate of disasters was about 3.09/year. It went down to about .93/year after the switchpoint.