Piecewise linear

[rstaubh]

I have the following signal which is piecewise linear and continuous. Currently I fail to use ruptures to extract my change points. The change points are known in advance. I have used clinear model and the Dynp.

Could you help me resolve the issue? Is there an alternative better solution to solve this?

# Define a piece-wise linear function
def piecewise_linear(x, x0, x1, m0, m1, m2):
    return np.piecewise(x, [x < x0, (x >= x0) & (x < x1), x >= x1], [lambda x: m0*x, lambda x: m1*(x - x0) + m0*x0, lambda x: m2*(x - x1) + m1*(x1 - x0) + m0*x0])
# Generate the true signal
n = 1000
x = np.arange(n)
x1 = 300
x2 = 700
m0 = 0.1
m1 = 0.4
m2 = 0.8
y = piecewise_linear(x, x1, x2, m0,m1,m2)

# Generate the noise
noise_level = 30
noise = noise_level * np.random.randn(n)

# Add noise to the signal
signal = y + noise
signal = signal.reshape((-1, 1))
# Plot the results
fig, ax = plt.subplots()
ax.plot(x, signal, label='Noisy signal')
ax.plot(x, y, label='True signal')
ax.legend()
plt.show()

# True breaking points:
bkps = [x1,x2,n]
algo = rpt.Dynp(model="clinear").fit(signal)
result = algo.predict(n_bkps=2)
rpt.display(signal, bkps, result)
print(f"{bkps=}")
print(f"{result=}")
plt.show()

[deepcharles]

Hi,

Currently, this cost funtion is not very robust. You can instead try to find changes in the slope without enforcing the continuity constraint. This has been shown to remain consistent. Basically, you replace the algo by

algo = rpt.Dynp(model="linear").fit(np.c_[signal.flatten(), x, np.ones_like(x)])

The full code is as follows:

# Define a piece-wise linear function
def piecewise_linear(x, x0, x1, m0, m1, m2):
    return np.piecewise(x, [x < x0, (x >= x0) & (x < x1), x >= x1], [lambda x: m0*x, lambda x: m1*(x - x0) + m0*x0, lambda x: m2*(x - x1) + m1*(x1 - x0) + m0*x0])
# Generate the true signal
n = 1000
x = np.arange(n)
x1 = 300
x2 = 700
m0 = 0.1
m1 = 0.4
m2 = 0.8
y = piecewise_linear(x, x1, x2, m0,m1,m2)

# Generate the noise
noise_level = 30
noise = noise_level * np.random.randn(n)

# Add noise to the signal
signal = y + noise
# Plot the results
fig, ax = plt.subplots()
ax.plot(x, signal, label='Noisy signal')
ax.plot(x, y, label='True signal')
ax.legend()
plt.show()

# True breaking points:
bkps = [x1,x2,n]
algo = rpt.Dynp(model="linear").fit(np.c_[signal.flatten(), x, np.ones_like(x)])
result = algo.predict(n_bkps=2)
fig, (ax,) = rpt.display(signal, bkps, result)
print(f"{bkps=}")
print(f"{result=}")
plt.show()

# reconstruction

from scipy.interpolate import splrep, BSpline
tck = splrep(x=x, y=signal.flatten(), k=1, t=result[:-1])
reconstruction = BSpline(*tck)(x)
fig, ax = plt.subplots()
ax.plot(x, signal, label='Noisy signal')
ax.plot(x, y, label='True signal')
ax.plot(x, reconstruction, label='reconstruction')
ax.legend()
plt.show()

[rstaubh]

Thanks for your answer!