Кусочно линейная подгонка с n точками останова

Я использовал код, найденный в вопросе Как применить кусочно-линейную посадку в Python? , для выполнения сегментированной линейной аппроксимации с единственной точкой останова.

Код выглядит следующим образом:

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)) 

Я пытаюсь выяснить, как я могу расширить это, чтобы обрабатывать n точек останова.

Я попробовал следующий код для метода кусочно_линейный () для обработки 2 точек останова, но он никак не изменяет значения точек останова.

 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)) с 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)) 

Любой вход был бы весьма полезен

    NumPy имеет функцию polyfit которая позволяет легко находить линию наилучшего соответствия через набор точек:

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

    Так что единственная трудность – найти точки останова. Для заданного набора точек останова тривиально найти наилучшие соответствия по данным.

    Поэтому вместо того, чтобы сразу найти местоположение точек останова и коэффициентов линейных частей, достаточно минимизировать пространство параметров точек останова.

    Так как точки останова могут быть заданы их целыми индексами в массив x , пространство параметров можно рассматривать как точки в целочисленной сетке из N измерений, где N – количество точек останова.

    optimize.curve_fit не является хорошим выбором в качестве минимизатора для этой проблемы, поскольку пространство параметров является целочисленным. Если вы будете использовать curve_fit , алгоритм будет настраивать параметры, чтобы определить, в каком направлении двигаться. Если настройка менее 1 единицы, значения x точек останова не изменяются, поэтому ошибка не изменяется, поэтому алгоритм не получает информации о правильном направлении, в котором необходимо изменить параметры. Следовательно, curve_fit имеет тенденцию терпеть неудачу, когда пространство параметров по существу целочисленно.

    Лучшим, но не очень быстрым, минимизатором будет поиск сетки грубой силы. Если число точек останова невелико (а пространство параметров x значений невелико), этого может быть достаточно. Если число точек останова велико и / или пространство параметров велико, то, возможно, настройте многостадийный поиск грубой / тонкой (грубой силы) сетки. Или, возможно, кто-то предложит более разумный минимизатор, чем грубая сила …


     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() 

    печать

     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] 

    и участки

    введите описание изображения здесь