import pymc as pm
import numpy as np
import arviz as az
from pymc.math import switch, ge, exp

%load_ext lab_black

3. Revisiting UK Coal Mining Disasters*#

Adapted from Unit 10: disasters.odc.

Data can be found here.

Change Point Analysis, discussed previously in this Unit 5 example about Gibbs sampling.

The 112 data points represent the numbers of coal-mining disasters involving 10 or more men killed per year between 1851 and 1962.

Based on the observation that the there was a significant decrease around 1900, it is suitable to apply a change-point model to divide the whole dataset into two periods; each period with its own distribution of number of disasters.

The data set was compiled by Maguire, Pearson and Wynn in 1952 and updated by Jarrett (1978). This data have been used by a number of authors to illustrate various techniques that can be applied to point processes

Maguire, B. A., Pearson, E. S. and Wynn, A. H. A. (1952). The time intervals between industrial accidents. Biometrika, 39, 168†180.

Jarrett, R.G. (1979). A note on the intervals between coal-mining disasters. Biometrika, 66, 191-193.

Carlin, Gelfand, and Smith (1992) Heirarchical Bayesian Analysis of Changepoint Problems. Applied Statistics, 41, 389-405.

# X is the number of coal mine disasters per year
# fmt: off
X = np.array([4, 5, 4, 1, 0, 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, 2, 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, 1, 2, 1, 1, 1, 1, 2, 4, 2,
     0, 0, 0, 1, 4, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])
# fmt: on

y = np.array([y for y in range(1851, 1963)])

Model 1#

α = 4
β = 1
γ = 0.5
δ = 1

with pm.Model() as m:
    year = pm.Uniform("year", 1851, 1963)
    λ = pm.Gamma("λ", α, β)
    μ = pm.Gamma("μ", γ, δ)

    diff = pm.Deterministic("diff", μ - λ)

    rate = λ + switch(ge(y - year, 0), 1, 0) * diff
    pm.Poisson("lik", mu=rate, observed=X)

    trace = pm.sample(2000)
Hide code cell output 
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [year, λ, μ]
100.00% [12000/12000 03:26<00:00 Sampling 4 chains, 0 divergences] 
Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 206 seconds.
Chain 0 reached the maximum tree depth. Increase `max_treedepth`, increase `target_accept` or reparameterize.
Chain 1 reached the maximum tree depth. Increase `max_treedepth`, increase `target_accept` or reparameterize.
Chain 2 reached the maximum tree depth. Increase `max_treedepth`, increase `target_accept` or reparameterize.
Chain 3 reached the maximum tree depth. Increase `max_treedepth`, increase `target_accept` or reparameterize.

az.summary(trace)
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
year 1890.405 2.405 1886.000 1894.461 0.108 0.076 537.0 500.0 1.0
λ 3.153 0.284 2.655 3.698 0.011 0.008 658.0 1025.0 1.0
μ 0.916 0.115 0.698 1.137 0.004 0.003 645.0 826.0 1.0
diff -2.237 0.301 -2.805 -1.681 0.012 0.008 670.0 1185.0 1.0

Model 2#

with pm.Model() as m:
    year = pm.Uniform("year", 1851, 1963)
    z0 = pm.Normal("z0", 0, tau=0.00001)
    z1 = pm.Normal("z1", 0, tau=0.00001)

    λ = pm.Deterministic("λ", exp(z0))
    μ = pm.Deterministic("μ", exp(z0 + z1))

    diff = pm.Deterministic("diff", μ - λ)

    rate = pm.math.exp(z0 + switch(ge(y - year, 0), 1, 0) * z1)
    pm.Poisson("lik", mu=rate, observed=X)

    trace = pm.sample(2000)
Hide code cell output 
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [year, z0, z1]
100.00% [12000/12000 05:01<00:00 Sampling 4 chains, 0 divergences] 
Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 302 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters.  A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details
Chain 0 reached the maximum tree depth. Increase `max_treedepth`, increase `target_accept` or reparameterize.
Chain 1 reached the maximum tree depth. Increase `max_treedepth`, increase `target_accept` or reparameterize.
Chain 2 reached the maximum tree depth. Increase `max_treedepth`, increase `target_accept` or reparameterize.
Chain 3 reached the maximum tree depth. Increase `max_treedepth`, increase `target_accept` or reparameterize.

az.summary(trace)
mean sd hdi_3% hdi_97% mcse_mean mcse_sd ess_bulk ess_tail r_hat
z0 1.143 0.089 0.971 1.306 0.005 0.003 378.0 748.0 1.01
z1 -1.233 0.145 -1.519 -0.972 0.006 0.004 529.0 1039.0 1.01
year 1890.258 2.311 1886.003 1893.954 0.117 0.083 409.0 478.0 1.01
λ 3.147 0.279 2.620 3.662 0.014 0.010 378.0 748.0 1.01
μ 0.921 0.116 0.703 1.141 0.004 0.003 946.0 1579.0 1.00
diff -2.227 0.289 -2.766 -1.689 0.015 0.011 368.0 683.0 1.01 
%load_ext watermark
%watermark -n -u -v -iv -p pytensor

Last updated: Wed Mar 22 2023

Python implementation: CPython
Python version       : 3.11.0
IPython version      : 8.9.0

pytensor: 2.10.1

pymc : 5.1.2
arviz: 0.15.1
numpy: 1.24.2