# Select optimal piecewise regression fit

I'm making a program which fits a piecewise linear regression with up to 4-5 breakpoints in the data, and then deciding how many breakpoints is best to prevent over and underfitting. However, my code is extremely slow to run, due to how ungraceful and brute force it is. However, I'm not sure how to code it in any other less laborious way then like so:

A rough draft I have of my code suiting my aim is this:

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

def segmentedRegression_two(xData,yData):

def func(xVals,break1,break2,slope1,offset1,slope_mid,offset_mid,slope2,offset2):
returnArray=[]
for x in xVals:
if x < break1:
returnArray.append(slope1 * x + offset1)
elif (np.logical_and(x >= break1,x<break2)):
returnArray.append(slope_mid * x + offset_mid)
else:
returnArray.append(slope2 * x + offset2)

return returnArray

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

return np.sum((yData-model_y)**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([-y_max,y_min])
initial_parameters.append([-slope,slope])
initial_parameters.append([-y_max,y_min])
initial_parameters.append([-slope,slope])
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)

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

return Rsquared

def segmentedRegression_three(xData,yData):

def func(xVals,break1,break2,break3,slope1,offset1,slope2,offset2,slope3,offset3,slope4,offset4):
returnArray=[]
for x in xVals:
if x < break1:
returnArray.append(slope1 * x + offset1)
elif (np.logical_and(x >= break1,x<break2)):
returnArray.append(slope2 * x + offset2)
elif (np.logical_and(x >= break2,x<break3)):
returnArray.append(slope3 * x + offset3)
else:
returnArray.append(slope4 * x + offset4)

return returnArray

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

return np.sum((yData-model_y)**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([x_max,x_min])
initial_parameters.append([-slope,slope])
initial_parameters.append([-y_max,y_min])
initial_parameters.append([-slope,slope])
initial_parameters.append([-y_max,y_min])
initial_parameters.append([-slope,slope])
initial_parameters.append([y_max,y_min])
initial_parameters.append([-slope,slope])
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)

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

return Rsquared

def segmentedRegression_four(xData,yData):

def func(xVals,break1,break2,break3,break4,slope1,offset1,slope2,offset2,slope3,offset3,slope4,offset4,slope5,offset5):
returnArray=[]
for x in xVals:
if x < break1:
returnArray.append(slope1 * x + offset1)
elif (np.logical_and(x >= break1,x<break2)):
returnArray.append(slope2 * x + offset2)
elif (np.logical_and(x >= break2,x<break3)):
returnArray.append(slope3 * x + offset3)
elif (np.logical_and(x >= break3,x<break4)):
returnArray.append(slope4 * x + offset4)
else:
returnArray.append(slope5 * x + offset5)

return returnArray

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

return np.sum((yData-model_y)**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([x_max,x_min])
initial_parameters.append([x_max,x_min])
initial_parameters.append([-slope,slope])
initial_parameters.append([-y_max,y_min])
initial_parameters.append([-slope,slope])
initial_parameters.append([-y_max,y_min])
initial_parameters.append([-slope,slope])
initial_parameters.append([y_max,y_min])
initial_parameters.append([-slope,slope])
initial_parameters.append([y_max,y_min])
initial_parameters.append([-slope,slope])
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)

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

return Rsquared

r2s=[segmentedRegression_two(xData,yData),segmentedRegression_three(xData,yData),segmentedRegression_four(xData,yData)]

best_fit=np.max(r2s)

best_fit

np.where(r2s==best_fit)

r2s

• @AlexV I don't see any issues with indentations for the code. Can you point out to me where the issues are? Jun 12, 2019 at 16:15
• Ah, I see what you mean. I'll look into this. Jun 12, 2019 at 16:16
• I've fixed it as well as I can to my knowledge. Jun 12, 2019 at 16:26
• The paper kinked below deals with this kind of problems. But the cases of "up to 4-5 breakpoints" isn't specifically treated. Nevertheless, the method shown in the paper might give you an hint. fr.scribd.com/document/380941024/… Sep 22, 2019 at 9:56

First of all you should definitely clean up your code. There are a lot of unused variables and duplicate code, likely caused by copy and paste. Since you did not provide test data, there is no chance to tell how much time would be gained simply by removing [R[M]]SE computation from the functions.

Apart from that, its likely that the different funcs are the main reason for the bad performance. In your current implementation your basically throwing everything away that allows NumPy to be fast: the ability to perform loops in C instead of in Python. For a more detailed introduction to this topic I would highly recommend to watch Losing your Loops Fast Numerical Computing with NumPy and Performance Python: Seven Strategies for Optimizing Your Numerical Code by Jake VanderPlas. There is also an online version of the Python Data Science Handbook also by Jake VanderPlas.

So how you could func (I chose the first one) be rewritten to allow NumPy to play out its strenghts? (Note: this code is untested since you did not provide test data)

def func(x_vals, break1, break2, slope1, offset1, slope_mid, offset_mid, slope2, offset2):
"""Regression fit function"""
y = np.zeros_like(x_vals)