Segmented regression in python using differential evolution

Question

I have the long-term aim of creating a module that for a specific data set, fits segmented regressions up to an arbitrary number of breakpoints, as well as a standard polynomial and linear curve fit, and then evaluates which of the fits are the most appropriate for the data (likely using AIC or BIC).

I have a function that uses differential evolution to use segmented regression on an x and y dataset assuming 1 breakpoint:

def segReg_one(xData,yData):

    def func(xVals,model_break,slopeA,slopeB,offsetA,offsetB): #Initialization of the piecewise function
        returnArray=[]
        for x in xVals:
            if x > model_break:
                returnArray.append(slopeA * x + offsetA)
            else:
                returnArray.append(slopeB * x + offsetB)


        return returnArray

    def sumSquaredError(parametersTuple): #Definition of an error function to minimize
        modely=func(xData,*parametersTuple)
        warnings.filterwarnings("ignore") # Ignore warnings by genetic algorithm

        return np.sum((yData-modely)**2.0)

    def generate_genetic_Parameters():
        initial_parameters=[]
        x_max=np.max(xData)
        x_min=np.min(xData)
        y_max=np.max(yData)
        y_min=np.min(yData)
        slope=10*(y_max-y_min)/(x_max-x_min)

        initial_parameters.append([x_max,x_min]) #Bounds for model break point
        initial_parameters.append([-slope,slope]) #Bounds for slopeA
        initial_parameters.append([-slope,slope]) #Bounds for slopeB
        initial_parameters.append([y_max,y_min]) #Bounds for offset A
        initial_parameters.append([y_max,y_min]) #Bounds for offset B

        result=differential_evolution(sumSquaredError,initial_parameters,seed=3)

        return result.x

    geneticParameters = generate_genetic_Parameters() #Generates genetic parameters



    fittedParameters, pcov= curve_fit(func, xData, yData, geneticParameters) #Fits the data 
    print('Parameters:', fittedParameters)
    print('Model break at: ', fittedParameters[0])
    print('Slope of line where x < model break: ', fittedParameters[1])
    print('Slope of line where x > model break: ', fittedParameters[2])
    print('Offset of line where x < model break: ', fittedParameters[3])
    print('Offset of line where x > model break: ', fittedParameters[4])





    model=func(xData,*fittedParameters)

    absError = model - yData

    SE = np.square(absError) 
    MSE = np.mean(SE) 
    RMSE = np.sqrt(MSE) 
    Rsquared = 1.0 - (np.var(absError) / np.var(yData))

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



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


            axes.plot(xData, yData,  'D')

            xModel = np.linspace(min(xData), max(xData))
            yModel = func(xModel, *fittedParameters)

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


    graphWidth = 800
    graphHeight = 600

    return ModelAndScatterPlot(800,600)

Which runs fine. However, I tried to expand the model to allow for more than 1 breakpoint:

def segReg_two(xData,yData):

    def func(xData,break1,break2,slope1,slope_mid,slope2,offset1,offset_mid,offset2):
        returnArray=[]
        for x in xData:
            if x < break1:
                returnArray.append(slope1 * x + offset1)

            if (x < break2 and x > break1):
                returnArray.append(slope_mid * x + offset_mid)

            else:
                returnArray.append(slope2 * x + offset2)


    def sumSquaredError(parametersTuple): #Definition of an error function to minimize
        modely=func(xData,*parametersTuple)
        warnings.filterwarnings("ignore") # Ignore warnings by genetic algorithm

        return np.sum((yData-modely)**2.0)

    def generate_genetic_Parameters():
        initial_parameters=[]
        x_max=np.max(xData)
        x_min=np.min(xData)
        y_max=np.max(yData)
        y_min=np.min(yData)
        slope=10*(y_max-y_min)/(x_max-x_min)

        initial_parameters.append([x_max,x_min]) #Bounds for model break point
        initial_parameters.append([x_max,x_min])
        initial_parameters.append([-slope,slope]) 
        initial_parameters.append([-slope,slope]) 
        initial_parameters.append([-slope,slope]) 
        initial_parameters.append([y_max,y_min])
        initial_parameters.append([y_max,y_min]) 
        initial_parameters.append([y_max,y_min]) 


        result=differential_evolution(sumSquaredError,initial_parameters,seed=3)

        return result.x

    geneticParameters = generate_genetic_Parameters() #Generates genetic parameters



    fittedParameters, pcov= curve_fit(func, xData, yData, geneticParameters) #Fits the data 
    print('Parameters:', fittedParameters)
    print('Model break at: ', fittedParameters[0])
    print('Slope of line where x < model break: ', fittedParameters[1])
    print('Slope of line where x > model break: ', fittedParameters[2])
    print('Offset of line where x < model break: ', fittedParameters[3])
    print('Offset of line where x > model break: ', fittedParameters[4])





    model=func(xData,*fittedParameters)

    absError = model - yData

    SE = np.square(absError) 
    MSE = np.mean(SE) 
    RMSE = np.sqrt(MSE) 
    Rsquared = 1.0 - (np.var(absError) / np.var(yData))

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

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


            axes.plot(xData, yData,  'D')

            xModel = np.linspace(min(xData), max(xData))
            yModel = func(xModel, *fittedParameters)

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


    graphWidth = 800
    graphHeight = 600

    return ModelAndScatterPlot(800,600)

And this code runs into problems when I run segReg_two(x,y), stopping at the differential_evolution bit:

TypeError: unsupported operand type(s) for -: 'float' and 'NoneType' During handling of the above exception, another exception occurred:

RuntimeError: The map-like callable must be of the form f(func, iterable), returning a sequence of numbers the same length as 'iterable'

I didn't have this problem with segReg_one, so I don't see why I'm having it happen here. I am assuming (and I may be incorrect with this assumption) that the argument iterable must have compatible dimensions with my error function. However, I'm not exactly sure on how those two arguments exactly relate other than the fact that I'm finding the breakpoints, slopes and offsets that minimize the breakpoints given the bounds I have.

Also, my plan of attack seems extremely long-winded and brutish. Is there a better way to tackle this?

I think perhaps it is considering my piecewise function as None-type. Printing the function with some random values returned simply "None". However, my piecewise function prints the same thing and it still worked out fine.

Segmented regression in python using differential evolution

Answers (1)

Related Questions