Change point detection in linear regression

[YungDurum]

Hey,

I just started using the ruptures module and I have a question related to this module.

My data consist of many lineair regressions. I use ruptures to detect the change points.

Because I'm interested in the slope of the lineair regressions , I use the 'pwlf' module to determine the slope.

For ruptures I use the following settings: search engine = Pelt, cost function = l1 (only one tested so far). The jump, penalty-value and min_size I vary.

Observing the linear regressions the search engine seems to detect change points with almost zero change in slope. Is there a way I can set the minimum change in slope to detect?

Thanks in advance,

[deepcharles]

the cost function l1 detects change in the mean (in the median actually), so it will certainly make mistakes in case of change of slopes.

You can try the cost function clinear.

Hope this helps.

[oboulant]

Hello @YungDurum ,

Indeed, as @deepcharles suggested, if in your data you have continuity at the slope change point, then clinear cost function might help you. Here is the link to the documentation for this cost function.

Hope it helps !

[YungDurum]

Hi oboulant, I can't open the link.

[oboulant]

Hi @YungDurum , sorry, my mistake. I corrected the link. The link should be okay now.

[oboulant]

@YungDurum, if you signal has discontinuity at the break points, you still might have some solutions. I see two possibilities depending on the shape of the signal you are trying the segment :

• If the shape of the signal you are trying to segment is 1-dimensional, you might be able to segment your signal using CostLinear with a simple np.arange(0, signal.shape[0]) as the X matrix. Beware of how to package the signal in order to pass both the y and the X at the cost object creation time : you can look at the example here and the CostLinear class doc to see how to stack both y and X
• If the shape of the signal you are trying to segment is N-dimensional, then you still can create your own custom cost class that inherits from the BaseCost class. In your own custom class, you can define in the error method the linear regression. Once you have your custom error class, then you can pass it to your chosen search method : see this for instance for an example on how to pass a custom error class to the search method

Hope it helps !

[oboulant]

If you know a priori the number of breakpoints 🔽

import matplotlib.pyplot as plt  # for display purposes
from numpy.random import RandomState
import numpy as np
import ruptures as rpt

# Create dummy piecewise linear signal with discontinuity at change points 
x = np.arange(1000)
signal = 1.0 + 1.0*x
signal2 = -2.0 + 8.0*x
signal3 = -4.0 - 2.0*x
signal4 = 2.0*np.ones((1000)) 
signal[200:400] = signal2[200:400]
signal[400:600] = signal3[400:600]
signal[600:] = signal4[600:]

# Add noise
rd = RandomState(1234567890)
signal = signal + rd.normal(scale=200.0, size=1000)

# Run ruptures algo
# 1. Choose Dynp to run the most accurate (and costly) algorythm
# 2. Choose the "linear" cost model
# 3. Stack `signal` with `x`, `CostLinear` needs it to run the linear regressions
algo = rpt.detection.Dynp(model="linear",min_size=1, jump=1).fit(np.column_stack((signal.reshape(-1, 1), x)))
bkps = algo.predict(n_bkps=3)
print(bkps)

# Display
figure, axs = rpt.display(signal, [200, 400, 600, 1000], bkps)
figure.show()

[oboulant]

If you do not know a priori the number of breakpoints 🔽

import matplotlib.pyplot as plt  # for display purposes
from numpy.random import RandomState
import numpy as np
import ruptures as rpt

# Create dummy piecewise linear signal with discontinuity at change points 
x = np.arange(1000)
signal = 1.0 + 1.0*x
signal2 = -2.0 + 8.0*x
signal3 = -4.0 - 2.0*x
signal4 = 2.0*np.ones((1000)) 
signal[200:400] = signal2[200:400]
signal[400:600] = signal3[400:600]
signal[600:] = signal4[600:]


# Add noise
rd = RandomState(1234567890)
signal = signal + rd.normal(scale=200.0, size=1000)

# Run ruptures algo
# 1. Choose Pelt to run the most accurate (and costly) algorythm
# 2. Choose the "linear" cost model
# 3. Stack `signal` with `x`, `CostLinear` needs it to run the linear regressions
# 4. Choose a penalty value 
algo = rpt.detection.Pelt(model="linear",min_size=1, jump=1).fit(np.column_stack((signal.reshape(-1, 1), x)))
bkps = algo.predict(pen=10**6)
print(bkps)

# Display
figure, axs = rpt.display(signal, [200, 400, 600, 1000], bkps)
figure.show()

[deepcharles]

Note that if you need faster (but slightly less accurate) results, you can set jump=5 (or more) to only consider indexes that are multiples of 5.

[YungDurum]

Thanks for helping me out! I noticed I misinterpret my data and my data is continuous. Will I still be able to use your example code? Also I noticed that the 'cost linear model' is sometimes referred to as "clinear" or "linear", do they refer to the same function?

[oboulant]

Hi,

The above code sample still works, but it will give you the best change points without taking care of the continuity constraints. The choice of the cost function really relies on the underlying assumption you make on your data.

If, given your data, the continuity at the change point is a structural constraints, here is a code example :

import numpy as np
import matplotlib.pylab as plt
import ruptures as rpt
from numpy.random import RandomState

# creation of data
signal = 2.0*np.ones((200)) 
signal2 = -3.0*np.ones((200)) 
signal3 = np.zeros((200)) 
signal4 = 1.0*np.ones((200)) 
signal = np.cumsum(np.concatenate((signal, signal2, signal3, signal4), axis = None))

# Add noise
rd = RandomState(123)
signal = signal + rd.normal(scale=20.0, size=800)

# Run ruptures algo
# 1. Choose Dynp to run the most accurate (and costly) algorythm
# 2. Choose the "clinear" cost model
algo = rpt.detection.Dynp(model="clinear",min_size=1, jump=1).fit(signal)
bkps = algo.predict(n_bkps=3)
print(bkps)

figure, axs = rpt.display(signal, [200, 400, 600, 800], bkps)
figure.show()

As for your comment 'cost linear model' is sometimes referred to as "clinear" or "linear", could you point to us where exactly ? What is for sure is that model="clinear" is different than model="linear".

Hope it helps !

[oboulant]

Hi @YungDurum ,

I do not know in which setup you are working on this, but just so you know if the problem to solve is research oriented, it is possible for us at Centre Borelli to work on a joint paper publication in a scientific journal.
Moreover, if your data is public, we would be happy to create an example based on it to be including in ruptures public documentation.

Enjoy using ruptures !

[YungDurum]

Hi Oboulant,

I'm currently working on my bachelor thesis for the Vrije Universiteit Amsterdam at the faculty Physics of Living systems. I will definitely refer in my group and to my supervisor about your institution. Once finished, I will send my paper and data analysis to you, and your free to use it.

Cheers!

[1a1a11a]

Thank you for all the examples! This is very useful!
Is there an easy way to retrieve the slope of each segment?

[deepcharles]

If you have a function that compute the slope, say compute_slope(), you could do

list_of_slopes = list()
for (start, end) in rpt.utils.pairwise([0] + bkps):
    slope = compute_slope(signal[start:end])
    list_of_slopes.append(slope)

where signal is the signal at hand and bkps is a list a change-point indexes.

[1a1a11a]

To make sure I understand, rupture does not provide the slope as an output even though its optimization uses the slope, is this correct understanding?

[deepcharles]

Yes indeed. It is computed but kept in memory.