Bibliothèque Python pour la régression segmentée (ou régression par morceaux)

16

Je recherche une bibliothèque Python qui peut effectuer une régression segmentée (ou régression par morceaux) .

Exemple :

entrez la description de l'image ici

Franck Dernoncourt
la source
Cette question donne une méthode pour effectuer une régression par morceaux en définissant une fonction et en utilisant des bibliothèques python standard. stackoverflow.com/questions/29382903/…
Une question similaire ( stackoverflow.com/questions/29382903/… ) et une bibliothèque utile pour la régression par morceaux ( pypi.org/project/pwlf )
prashanth

Réponses:

7

numpy.piecewise peut le faire.

par morceaux (x, condlist, funclist, * args, ** kw)

Évaluez une fonction définie par morceaux.

Étant donné un ensemble de conditions et de fonctions correspondantes, évaluez chaque fonction sur les données d'entrée partout où sa condition est vraie.

Un exemple est donné sur SO ici . Pour être complet, voici un exemple:

from scipy import optimize
import matplotlib.pyplot as plt
import numpy as np

x = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10 ,11, 12, 13, 14, 15], dtype=float)
y = np.array([5, 7, 9, 11, 13, 15, 28.92, 42.81, 56.7, 70.59, 84.47, 98.36, 112.25, 126.14, 140.03])

def piecewise_linear(x, x0, y0, k1, k2):
    return np.piecewise(x, [x < x0, x >= x0], [lambda x:k1*x + y0-k1*x0, lambda x:k2*x + y0-k2*x0])

p , e = optimize.curve_fit(piecewise_linear, x, y)
xd = np.linspace(0, 15, 100)
plt.plot(x, y, "o")
plt.plot(xd, piecewise_linear(xd, *p))
christopherlovell
la source
4

La méthode proposée par Vito MR Muggeo [1] est relativement simple et efficace. Il fonctionne pour un nombre spécifié de segments et pour une fonction continue. Les positions des points d'arrêt sont estimées de manière itérative en effectuant, pour chaque itération, une régression linéaire segmentée permettant des sauts aux points d'arrêt. A partir des valeurs des sauts, les positions de point d'arrêt suivantes sont déduites, jusqu'à ce qu'il n'y ait plus de discontinuité (sauts).

"le processus est itéré jusqu'à une éventuelle convergence, qui n'est pas, en général, garantie"

En particulier, la convergence ou le résultat peut dépendre de la première estimation des points d'arrêt.

Ceci est la méthode utilisée dans le R paquet segmenté .

Voici une implémentation en python:

import numpy as np
from numpy.linalg import lstsq

ramp = lambda u: np.maximum( u, 0 )
step = lambda u: ( u > 0 ).astype(float)

def SegmentedLinearReg( X, Y, breakpoints ):
    nIterationMax = 10

    breakpoints = np.sort( np.array(breakpoints) )

    dt = np.min( np.diff(X) )
    ones = np.ones_like(X)

    for i in range( nIterationMax ):
        # Linear regression:  solve A*p = Y
        Rk = [ramp( X - xk ) for xk in breakpoints ]
        Sk = [step( X - xk ) for xk in breakpoints ]
        A = np.array([ ones, X ] + Rk + Sk )
        p =  lstsq(A.transpose(), Y, rcond=None)[0] 

        # Parameters identification:
        a, b = p[0:2]
        ck = p[ 2:2+len(breakpoints) ]
        dk = p[ 2+len(breakpoints): ]

        # Estimation of the next break-points:
        newBreakpoints = breakpoints - dk/ck 

        # Stop condition
        if np.max(np.abs(newBreakpoints - breakpoints)) < dt/5:
            break

        breakpoints = newBreakpoints
    else:
        print( 'maximum iteration reached' )

    # Compute the final segmented fit:
    Xsolution = np.insert( np.append( breakpoints, max(X) ), 0, min(X) )
    ones =  np.ones_like(Xsolution) 
    Rk = [ c*ramp( Xsolution - x0 ) for x0, c in zip(breakpoints, ck) ]

    Ysolution = a*ones + b*Xsolution + np.sum( Rk, axis=0 )

    return Xsolution, Ysolution

Exemple:

import matplotlib.pyplot as plt

X = np.linspace( 0, 10, 27 )
Y = 0.2*X  - 0.3* ramp(X-2) + 0.3*ramp(X-6) + 0.05*np.random.randn(len(X))
plt.plot( X, Y, 'ok' );

initialBreakpoints = [1, 7]
plt.plot( *SegmentedLinearReg( X, Y, initialBreakpoints ), '-r' );
plt.xlabel('X'); plt.ylabel('Y');

graphique

[1]: Muggeo, VM (2003). Estimation des modèles de régression avec des points d'arrêt inconnus. Statistiques en médecine, 22 (19), 3055-3071.

xdze2
la source
3

Je cherchais la même chose, et malheureusement il semble qu'il n'y en ait pas pour le moment. Vous trouverez quelques suggestions sur la façon de procéder dans cette question précédente .

Alternativement, vous pouvez examiner certaines bibliothèques R, par exemple segmentées, SiZer, strucchange, et si quelque chose fonctionne, essayez d'incorporer le code R en python avec rpy2 .

Modification pour ajouter un lien vers py-earth , «Une implémentation Python des splines de régression adaptative multivariée de Jerome Friedman».

Saraw
la source
2

Il existe un article de blog avec une implémentation récursive de la régression par morceaux. Cette solution correspond à une régression discontinue.

Si vous n'êtes pas satisfait du modèle discontinu et que vous souhaitez un réglage continu, je vous propose de rechercher votre courbe dans une base de kcourbes en L, en utilisant Lasso pour la rareté:

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import Lasso
# generate data
np.random.seed(42)
x = np.sort(np.random.normal(size=100))
y_expected = 3 + 0.5 * x + 1.25 * x * (x>0)
y = y_expected + np.random.normal(size=x.size, scale=0.5)
# prepare a basis
k = 10
thresholds = np.percentile(x, np.linspace(0, 1, k+2)[1:-1]*100)
basis = np.hstack([x[:, np.newaxis],  np.maximum(0,  np.column_stack([x]*k)-thresholds)]) 
# fit a model
model = Lasso(0.03).fit(basis, y)
print(model.intercept_)
print(model.coef_.round(3))
plt.scatter(x, y)
plt.plot(x, y_expected, color = 'b')
plt.plot(x, model.predict(basis), color='k')
plt.legend(['true', 'predicted'])
plt.xlabel('x')
plt.ylabel('y')
plt.title('fitting segmented regression')
plt.show()

Ce code vous renverra un vecteur de coefficients estimés:

[ 0.57   0.     0.     0.     0.     0.825  0.     0.     0.     0.     0.   ]

En raison de l'approche Lasso, elle est clairsemée: le modèle a trouvé exactement un point de rupture parmi 10 possibles. Les nombres 0,57 et 0,825 correspondent à 0,5 et 1,25 dans le vrai DGP. Bien qu'elles ne soient pas très proches, les courbes ajustées sont:

entrez la description de l'image ici

Cette approche ne vous permet pas d'estimer exactement le point d'arrêt. Mais si votre jeu de données est suffisamment grand, vous pouvez jouer avec différents k(peut-être le régler par validation croisée) et estimer le point d'arrêt avec suffisamment de précision.

David Dale
la source