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.

  • 3
    $\begingroup$ Models that have a flexible number of parameters are often called "non-parametric". It's possible to build these in PyMC3, here is an example of how to do that: pymc-devs.github.io/pymc3/notebooks/dp_mix.html It's a bit more involved so I would recommend to start out simpler. $\endgroup$
    – twiecki
    Sep 12, 2016 at 8:40
  • $\begingroup$ Thank you, looks indeed way more complex than my current level of understanding :/ $\endgroup$
    – guest6012
    Sep 12, 2016 at 11:36
  • $\begingroup$ Do you actually need that model or just want to explore? In the former case you can just try a fixed number of switch-points. $\endgroup$
    – twiecki
    Sep 13, 2016 at 18:34
  • $\begingroup$ I am dealing in my day job with problems for which not specifying the number of 'switch-points' could potentially lead to interesting results, and thought it'll be useful to know. But yes, I can also see that it'll be easier to write code that evaluates several fixed scenarios, it'll just be less elegant and not as exhaustive I guess. $\endgroup$
    – guest6012
    Sep 14, 2016 at 11:48
  • $\begingroup$ Well, if you really want that, I'm sure it can be done using the non-parametric example I provided. $\endgroup$
    – twiecki
    Sep 15, 2016 at 13:53

1 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],
                   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])
  • $\begingroup$ Dear all, I have tried smba example with some random data of mine (that is Beta distributed), and an error message appeared: "ValueError: Unused step method arguments: {'steps'}". So, what would be t he cause of this message, and how to solve it? $\endgroup$ Nov 8, 2021 at 17:50
  • $\begingroup$ cheers, pymc3 has seen a lot of development during the past two years, without further context or an example, it's difficult to give advice here. $\endgroup$
    – smba
    Nov 11, 2021 at 14:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.