Ajustement linéaire par morceaux avec n points d'arrêt

Question

J'ai utilisé un code trouvé dans la question How to apply piecewise linear fit in Python?, pour effectuer une approximation linéaire segmentée avec un seul point d'arrêt.

Le code est le suivant:

from scipy import optimize 
import matplotlib.pyplot as plt 
import numpy as np 
%matplotlib inline 

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ,11, 12, 13, 14, 15], dtype=float) 
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03]) 

def piecewise_linear(x, x0, y0, k1, k2): 
    return np.piecewise(x, 
         [x < x0], 
         [lambda x:k1*x + y0-k1*x0, lambda x:k2*x + y0-k2*x0]) 

p , e = optimize.curve_fit(piecewise_linear, x, y) 
xd = np.linspace(0, 15, 100) 
plt.plot(x, y, "o") 
plt.plot(xd, piecewise_linear(xd, *p)) 

J'essaie de comprendre comment je peux l'étendre à gérer n points d'arrêt.

J'ai essayé le code suivant pour la méthode piecewise_linear() pour gérer 2 points d'arrêt, mais cela ne modifie en rien les valeurs des points d'arrêt.

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20], dtype=float) 
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03, 150, 152, 154, 156, 158]) 

def piecewise_linear(x, x0, x1, a1, b1, a2, b2, a3, b3): 
    return np.piecewise(x, 
         [x < x0, np.logical_and(x >= x0, x < x1), x >= x1 ], 
         [lambda x:a1*x + b1, lambda x:a2*x+b2, lambda x: a3*x + b3]) 

p , e = optimize.curve_fit(piecewise_linear, x, y) 
xd = np.linspace(0, 20, 100) 
plt.plot(x, y, "o") 
plt.plot(xd, piecewise_linear(xd, *p)) 

Toute entrée serait grandement apprécié

Answer 1

NumPy a une polyfit function ce qui le rend très facile de trouver la meilleure ligne d'ajustement par un ensemble de points:

coefs = npoly.polyfit(xi, yi, 1) 

Alors, vraiment la seule difficulté est de trouver les points d'arrêt. Pour un ensemble donné de points d'arrêt , il est trivial de trouver les meilleures lignes d'ajustement à travers les données données.

Ainsi, au lieu d'essayer de trouver l'emplacement des points d'arrêt et les coefficients des parties linéaires à la fois, il suffit de réduire au minimum sur un espace de paramètre de points d'arrêt.

Etant donné que les points d'arrêt peuvent être spécifiés par leurs valeurs d'index d'entier dans la gamme x, l'espace des paramètres peut être considéré comme des points sur une grille de nombre entier de N dimensions, où N est le nombre de points d'arrêt.

optimize.curve_fit n'est pas un bon choix en tant que minimiseur pour ce problème parce que l'espace de paramètre est de valeur entière. Si vous deviez utiliser curve_fit, l'algorithme modifierait les paramètres pour déterminer dans quelle direction déplacer. Si le tweak est inférieur à 1 unité, les valeurs x des points d'arrêt ne changent pas , donc l'erreur ne change pas, donc l'algorithme ne reçoit aucune information sur la direction correcte dans laquelle déplacer les paramètres. Par conséquent curve_fit tend à échouer lorsque l'espace de paramètre est essentiellement de valeur entière.

Un meilleur, mais pas très rapide, minimiseur serait une recherche de grille brute-force. Si le nombre de points d'arrêt est petit (et l'espace de paramètre de x -valeurs est petit) cela pourrait suffire. Si le nombre de points d'arrêt est important et/ou l'espace de paramètre est grand, alors peut-être mettre en place une recherche de grille grossière/fine multi-étape (force brute). Ou, peut-être que quelqu'un proposera un minimiseur plus intelligent que la force brute ...

---

import numpy as np 
import numpy.polynomial.polynomial as npoly 
from scipy import optimize 
import matplotlib.pyplot as plt 
np.random.seed(2017) 

def f(breakpoints, x, y, fcache): 
    breakpoints = tuple(map(int, sorted(breakpoints))) 
    if breakpoints not in fcache: 
     total_error = 0 
     for f, xi, yi in find_best_piecewise_polynomial(breakpoints, x, y): 
      total_error += ((f(xi) - yi)**2).sum() 
     fcache[breakpoints] = total_error 
    # print('{} --> {}'.format(breakpoints, fcache[breakpoints])) 
    return fcache[breakpoints] 

def find_best_piecewise_polynomial(breakpoints, x, y): 
    breakpoints = tuple(map(int, sorted(breakpoints))) 
    xs = np.split(x, breakpoints) 
    ys = np.split(y, breakpoints) 
    result = [] 
    for xi, yi in zip(xs, ys): 
     if len(xi) < 2: continue 
     coefs = npoly.polyfit(xi, yi, 1) 
     f = npoly.Polynomial(coefs) 
     result.append([f, xi, yi]) 
    return result 

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 
       18, 19, 20], dtype=float) 
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 
       126.14, 140.03, 150, 152, 154, 156, 158]) 
# Add some noise to make it exciting :) 
y += np.random.random(len(y))*10 

num_breakpoints = 2 
breakpoints = optimize.brute(
    f, [slice(1, len(x), 1)]*num_breakpoints, args=(x, y, {}), finish=None) 

plt.scatter(x, y, c='blue', s=50) 
for f, xi, yi in find_best_piecewise_polynomial(breakpoints, x, y): 
    x_interval = np.array([xi.min(), xi.max()]) 
    print('y = {:35s}, if x in [{}, {}]'.format(str(f), *x_interval)) 
    plt.plot(x_interval, f(x_interval), 'ro-') 


plt.show() 

impressions

y = poly([ 4.58801083 2.94476604]) , if x in [1.0, 6.0] 
y = poly([-70.36472935 14.37305793]) , if x in [7.0, 15.0] 
y = poly([ 123.24565235 1.94982153]), if x in [16.0, 20.0] 

et parcelles