Автоматическое определение кусочно-регрессионных точек останова в данных

Question

Мне поручено, например, получить следующие данные:

(оси не важны, а также то, что представляют данные для моеговопрос).Предполагая, что мои степени свободы таковы, что при подгонке кривых я не ожидаю смоделировать что-либо выше, чем с полиномами степени 3, и, как правило, все они должны быть линейными.

Правильное соответствие данных, по-видимому, выглядит следующим образом:

Однако мне удалось уместить эти две строки только потому, что япроизвольно отметил, что, как представляется, произошел разрыв в качественном поведении данных (переход от приблизительно постоянной к приблизительно линейной) в точке, указанной на моем графике, где линия изменила свой наклон, так как модели были выбраны произвольно линейными

Есть ли способ, которым это можно сделать без априорного знания?Я пытаюсь закодировать программу, которая автоматически открывает файл данных и выбирает точку между двумя кривыми.Я ищу руководства или объяснения или что-то подобное, что дало бы мне знания о том, как найти точку между двумя различными кривыми подгонками или моделями, а затем выбрал подходящую модель (с соответствующей степенью полинома) для каждого раздела данных..

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

Приложение

Ответ Джеймса Филипса выглядит как раз то, что я ищу, однако то, что яТеперь нужно понять, как работает код.Сейчас я отмечу мои замешательства с ним, так как ответы на эти вопросы теперь будут зависеть от того, как понять работу кода.

def sumOfSquaredError(parameterTuple):
    warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
    val = func(xData, *parameterTuple)
    return numpy.sum((yData - val) ** 2.0)

1. Что именно является val?Почему мы используем сумму квадратов ошибок, а не среднее значение, например?

def generate_Initial_Parameters(): # min and max used for bounds

    maxX = max(xData)
    minX = min(xData)
    maxY = max(yData)
    minY = min(yData)
    slope = 10.0 * (maxY - minY) / (maxX - minX) # times 10 for safety margin

    parameterBounds = []
    parameterBounds.append([minX, maxX]) # search bounds for breakpoint
    parameterBounds.append([-slope, slope]) # search bounds for slopeA
    parameterBounds.append([minY, maxY]) # search bounds for offsetA
    parameterBounds.append([-slope, slope]) # search bounds for slopeB
    parameterBounds.append([minY, maxY]) # search bounds for offsetB


    result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
    return result.x

Разве это не просто запрос данных x и yи затем создание массива с именем parameterBounds, который просто добавляет к нему массивы 1x2, кодирующие минимальные и максимальные значения x и y соответственно и два массива с наклоном в них.

1. Почему [-slope,slope] массивы добавляются к «границам поиска для slopeA» и «границам поиска для slopeB»?Почему он построен с использованием максимальных и минимальных значений x и y?Угадают ли они, чтобы найти правильные уклоны каждой части?

2. Что там делают 10 и почему это «запас прочности»?

3. Разве differential_evolution не используется для поиска минимума многомерной функции?Что именно он делает?

   modelPredictions = func(xData, *fittedParameters)

   absError = modelPredictions - yData

4. Как получилось, xData и *fittedParameters выполнили необходимое количествоаргументы для func?

Наконец, как этот код находит точку останова?

Answer 1

Вот пример кода, подгоняющего две разные прямые линии к набору данных, а также автоматически подгоняющего точку останова между двумя линиями. В этом примере используется стандартный модуль генетического алгоритма scipy absolute_evolution, в котором используется алгоритм Latin Hypercube для обеспечения тщательного поиска в пространстве параметров, требующего границ для поиска. В этом примере эти границы взяты из данных max и min.

import numpy, scipy, matplotlib
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
from scipy.optimize import differential_evolution
import warnings

xData = numpy.array([19.1647, 18.0189, 16.9550, 15.7683, 14.7044, 13.6269, 12.6040, 11.4309, 10.2987, 9.23465, 8.18440, 7.89789, 7.62498, 7.36571, 7.01106, 6.71094, 6.46548, 6.27436, 6.16543, 6.05569, 5.91904, 5.78247, 5.53661, 4.85425, 4.29468, 3.74888, 3.16206, 2.58882, 1.93371, 1.52426, 1.14211, 0.719035, 0.377708, 0.0226971, -0.223181, -0.537231, -0.878491, -1.27484, -1.45266, -1.57583, -1.61717])
yData = numpy.array([0.644557, 0.641059, 0.637555, 0.634059, 0.634135, 0.631825, 0.631899, 0.627209, 0.622516, 0.617818, 0.616103, 0.613736, 0.610175, 0.606613, 0.605445, 0.603676, 0.604887, 0.600127, 0.604909, 0.588207, 0.581056, 0.576292, 0.566761, 0.555472, 0.545367, 0.538842, 0.529336, 0.518635, 0.506747, 0.499018, 0.491885, 0.484754, 0.475230, 0.464514, 0.454387, 0.444861, 0.437128, 0.415076, 0.401363, 0.390034, 0.378698])


def func(xArray, breakpoint, slopeA, offsetA, slopeB, offsetB):
    returnArray = []
    for x in xArray:
        if x < breakpoint:
            returnArray.append(slopeA * x + offsetA)
        else:
            returnArray.append(slopeB * x + offsetB)
    return returnArray


# function for genetic algorithm to minimize (sum of squared error)
def sumOfSquaredError(parameterTuple):
    warnings.filterwarnings("ignore") # do not print warnings by genetic algorithm
    val = func(xData, *parameterTuple)
    return numpy.sum((yData - val) ** 2.0)


def generate_Initial_Parameters():
    # min and max used for bounds
    maxX = max(xData)
    minX = min(xData)
    maxY = max(yData)
    minY = min(yData)
    slope = 10.0 * (maxY - minY) / (maxX - minX) # times 10 for safety margin

    parameterBounds = []
    parameterBounds.append([minX, maxX]) # search bounds for breakpoint
    parameterBounds.append([-slope, slope]) # search bounds for slopeA
    parameterBounds.append([minY, maxY]) # search bounds for offsetA
    parameterBounds.append([-slope, slope]) # search bounds for slopeB
    parameterBounds.append([minY, maxY]) # search bounds for offsetB


    result = differential_evolution(sumOfSquaredError, parameterBounds, seed=3)
    return result.x

# by default, differential_evolution completes by calling curve_fit() using parameter bounds
geneticParameters = generate_Initial_Parameters()

# call curve_fit without passing bounds from genetic algorithm
fittedParameters, pcov = curve_fit(func, xData, yData, geneticParameters)
print('Parameters:', fittedParameters)
print()

modelPredictions = func(xData, *fittedParameters) 

absError = modelPredictions - yData

SE = numpy.square(absError) # squared errors
MSE = numpy.mean(SE) # mean squared errors
RMSE = numpy.sqrt(MSE) # Root Mean Squared Error, RMSE
Rsquared = 1.0 - (numpy.var(absError) / numpy.var(yData))

print()
print('RMSE:', RMSE)
print('R-squared:', Rsquared)

print()


##########################################################
# graphics output section
def ModelAndScatterPlot(graphWidth, graphHeight):
    f = plt.figure(figsize=(graphWidth/100.0, graphHeight/100.0), dpi=100)
    axes = f.add_subplot(111)

    # first the raw data as a scatter plot
    axes.plot(xData, yData,  'D')

    # create data for the fitted equation plot
    xModel = numpy.linspace(min(xData), max(xData))
    yModel = func(xModel, *fittedParameters)

    # now the model as a line plot
    axes.plot(xModel, yModel)

    axes.set_xlabel('X Data') # X axis data label
    axes.set_ylabel('Y Data') # Y axis data label

    plt.show()
    plt.close('all') # clean up after using pyplot

graphWidth = 800
graphHeight = 600
ModelAndScatterPlot(graphWidth, graphHeight)