Слишком медленный код при выборе оптимального кусочно-регрессионного соответствия

Question

Я делаю программу, которая подгоняет кусочно-линейную регрессию к 4-5 точкам останова в данных, а затем решаю, сколько точек останова лучше всего предотвратить, чтобы избежать чрезмерного или недостаточного соответствия. Тем не менее, мой код крайне медленен для запуска из-за того, насколько он изящен.

Черновик моего кода:

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)

Хотя мне может понадобиться AIC или что-то в этом роде.

Есть ли способ, как я могу сделать это более эффективным при работе?

Answer 1

Я взял один из ваших func и поместил его в тестовый скрипт:

import numpy as np

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

arr = np.linspace(0,20,10000)
breaks = [4, 10, 15]
slopes = [.1, .2, .3, .4]
offsets = [1,2,3,4]
sl_off = np.array([slopes,offsets]).T.ravel().tolist()
print(sl_off)
ret = func(arr, *breaks, *sl_off)
if len(ret)<25:
    print(ret)

Затем я сделал первый шаг «векторизации», оценивая функцию в блоках значений, а не элемент за элементом.

def func1(xVals, breaks, slopes, offsets):
    res = np.zeros(xVals.shape)
    i = 0 
    mask = xVals<breaks[i]
    res[mask] = slopes[i]*xVals[mask]+offsets[i]
    for i in [1,2]:
        mask = np.logical_and(xVals>=breaks[i-1], xVals<breaks[i])
        res[mask] = slopes[i]*xVals[mask]+offsets[i]
    i=3
    mask = xVals>=breaks[i-1]
    res[mask] = slopes[i]*xVals[mask]+offsets[i]
    return res

ret1 = func1(arr, breaks, slopes, offsets)
print(np.allclose(ret, ret1))

allclose тестовые распечатки True. Я также ran это в ipython, и приурочен к двум версиям.

In [41]: timeit func(arr, *breaks, *sl_off)                                                            
66.2 ms ± 337 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
In [42]: timeit func1(arr, breaks, slopes, offsets)                                                    
165 µs ± 586 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Я также сделал plt.plot(xVals, ret), чтобы увидеть простой график функции.

Я написал func1 с целью заставить его работать для всех 3 ваших дел. Его нет, но его не сложно изменить, чтобы оно зависело от длины входных списков (или массивов).

Я уверен, что можно сделать больше, но это должно начаться в правильном направлении.

Существует также numpy piecewise оценщик:

np.piecewise(x, condlist, funclist, *args, **kw)

но мне кажется, что создание двух входных списков будет такой же работой.