How to build a PyMC model to detect multiple 'switch points'?

Question

In 'Bayesian Methods for Hackers' first chapter, Cam Davidson-Pilon presents an example model for detecting at what time point did a user's frequency of text-messaging changed.

This model assumes a single time point of change. How would you model a more generic solution where the number of change points (and their location) is to be found? (the result would be a probability distribution for the number of change points and distributions for each of the points' time)

My thought (total new to this) was to update the 'lambda_' function by adding a random variable for the number of points, and then based on the variable's value to slice the 'out' array and generate matching lambda variables. However I am not sure how to do that in terms of creating dynamic random variables and whether such a thing will collide with the function's 'deterministic' declaration.

Answer 1

It's been three years, but I believe this might be the approach mentioned in the comments by @twiecki. It uses a truncated Dirichlet Process Mixture Model to detect multiple change points without any prior knowledge of their number.

I have tried to disassemble their implementation and re-implement it on my own, so take the following with a grain of salt. Let's say we have a time-discrete time-series and we wish to extract change points, in particular the (most likely) number of change points and their location. For now, we are interested in change points in the mean of the time-series rather than changes in the variance, so let's assume the variance remains constant across the entire time series.

We start with a time-series with three segments (means: 1, 2, and 3) and two change points at locations 50, 100, respectively.

import numpy as np
import pymc3 as pm

signal = np.concatenate([
   np.full((50,), 1.0), 
   np.full((50,), 2.0), 
   np.full((50,), 3.0)
])

First, we transform our time-series into something that describes change in the time-series, so we compute something like the gradient for each point.

growth = np.abs(np.gradient(signal))
growth = growth / np.sum(growth) #normalize

From here, we can conceive growth as a measure of how much change is attributed to each point in time. Or, more mathematically, as the empirical density function of the probability that a time point is actually a change point. In our small example, the density function has two modes at 50 and 100.

You can now try to model the transformed data with a mixture model. The simplest approach would be to test mixture models with 1...k components and figure out the best fit with something like the Elbow method. Since we do not know/want to specify the number of change points (or components, in this context), you can resort to truncated mixture models as described in the comments above. You might want to start with something like this:

K = 15 # at most 15 components
N = growth.shape[0]

with pm.Model():

   # mixture model weights
   w = pm.Dirichlet('w', np.ones(K))

   component = pm.Categorical('component', w, shape=N)

   # independent means and precision for each component
   mu = pm.Uniform('mu', 0., 1.0 * N, shape=K)
   tau = pm.HalfNormal('tau', 1.0, shape=K)

   mix = pm.Normal('obs', mu=mu[component],
                   tau=tau[component],
                   shape=N, observed=growth)

   step1 = pm.Metropolis(vars=[w, mu, tau, mix])
   step2 = pm.ElemwiseCategorical([component], np.arange(K))

   trace = pm.sample(5000, steps=[step1, step2])