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])
