In [1]:
from pymc3 import *

import theano.tensor as t
from numpy import arange, array, ones, concatenate
from numpy.random import randint

__all__ = ['disasters_data', 'switchpoint', 'early_mean', 'late_mean', 'rate',
'disasters']

# Time series of recorded coal mining disasters in the UK from 1851 to 1962
disasters_data = array([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, 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, 1, 4,
0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1])
years = len(disasters_data)

In [2]:
with Model() as model:

# Prior for distribution of switchpoint location
switchpoint = DiscreteUniform('switchpoint', lower=0, upper=years)
# Priors for pre- and post-switch mean number of disasters
early_mean = Exponential('early_mean', lam=1.)
late_mean = Exponential('late_mean', lam=1.)

# Allocate appropriate Poisson rates to years before and after current
# switchpoint location
idx = arange(years)
rate = switch(switchpoint >= idx, early_mean, late_mean)

# Data likelihood
disasters = Poisson('disasters', rate, observed=disasters_data)

Applied log-transform to early_mean and added transformed early_mean_log to model.
Applied log-transform to late_mean and added transformed late_mean_log to model.

In [3]:
%matplotlib inline

with model:

# Initial values for stochastic nodes
start = {'early_mean': 2., 'late_mean': 3.}

tr = sample(1000, start=start)

traceplot(tr, varnames=[early_mean, late_mean], priors=[Exponential.dist(1)]*2)

Assigned Metropolis to switchpoint
Assigned NUTS to early_mean_log
Assigned NUTS to late_mean_log
[-----------------100%-----------------] 1000 of 1000 complete in 0.8 sec