2
$\begingroup$

I have a series (sensor readings) with different events and I'm trying to segment/detect those events for further analysis using change point detection in python.

In general, the series consists of a mean with a fairly small standard deviation. Occasionally, an outlier reading will occur for one, maybe two points, and then usually revert to the same mean. Another common pattern is the series shifting to a new mean, either abruptly or over the course of several readings. I'd like the method I use to capture both kinds of changes.

Additionally, I'm going to be doing this online so I'd like the results to be somewhat stable. It's fine if new data makes a mean-shift more obvious and a new change point is detected at N-3, but I don't want a method that causes a change point to be detected at N-20.

It seems like ruptures is the most common change point library in python (I also looked at changefinder), all the methods are offline but I'd be happy re-running for every new data point if it worked. However, I've had trouble getting satisfactory results.

import ruptures as rpt
import matplotlib.pyplot as plt

signal = [9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 28.0, 13.0, 9.0, 10.0, 10.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 10.0, 10.0, 10.0, 9.0, 9.0, 9.0, 9.0, 9.0, 10.0, 31.0, 31.0, 35.0, 35.0, 37.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0]

def detect_events(signal):
    signal = np.array(signal)
    algo = rpt.Pelt(model="rbf", min_size=1, jump=1).fit(signal)
    result = algo.predict(pen=1)

    algo = rpt.Window(model='l2', width=2, min_size=1).fit(signal)
    result = algo.predict(pen=5)

    # algo = rpt.Binseg(model="l2", jump=10).fit(signal)
    # result = algo.predict(pen=20)

    change_points = [i for i in result if i < len(signal)]
    return change_points

change_points = detect_events(signal)

optimal_change_points = [7,9, 38, 43]

print(change_points)
for cp in change_points:
    print(signal[cp-2:cp+2])
f, ax = plt.subplots(1,1)
plt.plot(signal, axes=ax)

for cp in change_points:
    ax.axvline(x=cp, color='r', linestyle='--')

for cp in optimal_change_points:
    ax.axvline(x=cp, color='g', linestyle=':')

Series with optimal change points vs detected changepoints

It seems whatever combination of search method I try with ruptures I'm either missing significant events, or detecting loads of spurious events, or getting a change point at weird spots, like well into the sequence of -1 instead of at the start.

I'm getting to the point where I'm probably going to make my own method (I'm primarily concerned with mean-shifts so I think it's feasible), though I can see that turning into a giant pile of edge-cases. The problem feels standard enough that I'd just be re-inventing a much more functional wheel but I don't know where that wheel is.

$\endgroup$

1 Answer 1

1
$\begingroup$

Hello again (I commented on your question on Maths SE). I think you can get more desirable results using a different cost function for the set of change points rather than completely starting again with your own method.

Here is my own coding of PELT in Python. To start, I can recover something like the segmentation that you're describing if I run

import PELT

x = [9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 28.0, 13.0, 9.0, 10.0, 10.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 9.0, 10.0, 10.0, 10.0, 9.0, 9.0, 9.0, 9.0, 9.0, 10.0, 31.0, 31.0, 35.0, 35.0, 37.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0, -1.0]

PELT.PELT_cp_n(x, len(x) - 1, penalty = "BIC")[-1]

PELT_cp_n outputs the optimal set of change points that was determined at each time point, so [-1] selects the final set of change points. I get output

[0, 7, 8, 12, 29, 32, 37, 38, 40, 42, 43]

which looks like this: too many change points I agree with you that this is too many change points, particularly the segmentation of the long sequence of 9s and 10s, and it's happening because PELT can get a perfect fit on the long sequence of 9s, and it's actually more expensive for it to include a few 10s in the segment than to start a new segment for the 10s. Because I'm using the BIC penalty, every time PELT throws down a new change point it has to pay the complexity penalty $\beta = 2\log(n)$, where $n$ is the number of data points. If you increase the penalty $\beta$ to be more strict then you can enforce a sparser model that only triggers more serious change points. For example, if you change the start of the PELT_cp_n to be

def PELT_cp_n(x, end_time, likelihood_model = "gaussian", K = 0, penalty = "BIC"):
    if penalty == "BIC":
        beta = 2.0 * np.log(end_time + 1)
    elif penalty == "AIC":
        beta = 2.0 * 2.0
    elif penalty == "aluchko":
        beta = 25

and run

PELT.PELT_cp_n(x, len(x) - 1, penalty = "aluchko")

then you get output [0, 7, 9, 38, 43], depicted below. You could experiment with different penalties until you find one that works in general for the kind of data that you have.

right number of change points

$\endgroup$
1
  • 1
    $\begingroup$ Thanks for the help, that definitely gets me a lot closer. A constant penalty makes sense since it shouldn't be based on N. But another thought is I want to make it somewhat resistant to a shift in mean. For instance, a 20 surrounded by 19s x = [19] * 100 + [20] + [19] * 50 That's not a significant change for the phenomena but is a huge outlier for the Gaussian. I modified your implementation to add an offset to the standard deviation and that seems to perform better sigma_squared_hat = np.std(self.x[t0:t1]) + 1 It feels very hacky but it seems to work. $\endgroup$
    – aluchko
    Commented May 31, 2023 at 4:55

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.