Оптимизация кусочно-линейной регрессии

Question

Я написал функцию, которая при заданных параметрах может применять кусочно-линейную аппроксимацию с произвольным количеством кусочных секций к некоторым данным.

Я пытаюсь подогнать функцию под свои данные с помощью scipy.optimize.curve_fit, но получаю ошибку «OptimizeWarning: Ковариация параметров не может быть оценена». Я полагаю, что это может быть из-за вложенных лямбда-функций, которые я использую для определения кусочных секций.

Есть ли простой способ настройки моего кода, чтобы обойти это, или другая функция оптимизации scipy, которая может быть более подходящей?

#The piecewise function
def piecewise_linear(x, *params):

    N=len(params)/2

    if N.is_integer():N=int(N)
    else:raise(ValueError())

    c=params[0]
    xbounds=params[1:N]
    grads=params[N:]


    #First we define our conditions, which are true if x is a member of a given
    #bin.
    conditions=[]
    #first and last bins are a special case:
    cond0=lambda x: x<xbounds[0]
    condl=lambda x: x>=xbounds[-1]
    conditions.append(cond0(x))

    for i in range(len(xbounds)-1):
        cond=lambda x : (x >= xbounds[i]) & (x < xbounds[i+1])
        conditions.append(cond(x))

    conditions.append(condl(x))

    #Next we define our linear regression function for each bin. The offset
    #for each bin depends on where the previous bin ends, so we define
    #the regression functions recursively:

    functions=[]
    func0 = lambda x: grads[0]*x +c

    functions.append(func0)

    for i in range(len(grads)-1):
        func = (lambda j: lambda x: grads[j+1]*(x-xbounds[j])\
               +functions[j](xbounds[j]))(i)

        functions.append(func)

    return np.piecewise(x,conditions,functions)
#Some data

x=np.arange(100)
y=np.array([*np.arange(0,19,1),*np.arange(20,59,2),\
*np.arange(60,20,-1),*np.arange(21,42,1)]) + np.random.randn(100)

#A first guess of parameters
cguess=0
boundguess=[20,30,50]
gradguess=[1,1,1,1]
p0=[cguess,*boundguess,*gradguess]

fit=scipy.optimize.curve_fit(piecewise_linear,x,y,p0=p0)

Answer 1

Вот пример кода, который вписывает две прямые линии в криволинейный набор данных с точкой останова, где все параметры линии и точка останова установлены.В этом примере используется генетический алгоритм Scipy's diffrential Evolution для определения начальных оценок параметров регрессии.Этот модуль использует алгоритм Латинского гиперкуба, чтобы обеспечить тщательный поиск пространства параметров, для которого требуются границы, в которых производится поиск.В этом примере эти границы поиска получены из самих данных.Обратите внимание, что гораздо проще найти диапазоны для начальных оценок параметров, чем давать конкретные значения.

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)