分段线性函数的拐点检测

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

N = 100
s = [10 * i for i in range(1,10)]
delta = [0.1, 0,0,0,-5,0,0,11,5]
km = np.array([1,10])

t = np.linspace(0, 1, N)
t1 = np.ones((N, 2))
t1[:,0] = t

A = np.zeros((N, len(s)))
for i in range(len(s)):
    A[s[i]:,i] = t[:-s[i]]

y = t1 @ km +  A @ np.array(delta)

y_true = y
plt.plot(t, y_true )

def find_changepoint(y):
    N = len(y)
    NCP = 20
    s = (np.linspace(0,1,NCP+2)[1:-1]*N).astype(int)
    t = np.linspace(0, 1, N)
    t1 = np.ones((N, 2))
    t1[:,0] = t
    A = np.zeros((N, len(s)))           
    for i in range(len(s)):
        A[s[i]:,i] = t[:-s[i]]
    A = np.hstack([t1, A])

    from functools import partial
    from scipy.optimize import minimize
    
    def loss(A, y, delta):
        return np.sum(np.square(y - A @ delta)) + 0.001*np.sum(np.abs(delta))

    f = partial(loss, A, y)

    d0 = np.random.randn(len(s)+2)
    res = minimize(f, x0=d0)
    delta = res.x
        
    return delta, A @ delta

        
d, y_pred = find_changepoint(y_true)
plt.subplot(1,2,1)
plt.bar(range(len(delta)),delta, alpha=0.5)
plt.subplot(1,2,2)
plt.bar(range(len(d)-2),d[2:], alpha=0.5)
plt.figure()
plt.plot(t, y_pred, label='fit')
plt.plot(t, y_true, label='true')
plt.legend()

如下左图为真实的4个拐点的斜率变化，右图为检测到的斜率变化