Pustaka Python untuk regresi tersegmentasi (alias regresi satu demi satu)

Saya mencari pustaka Python yang dapat melakukan regresi tersegmentasi (alias regresi satu demi satu ) .

Contoh :

numpy.piecewise dapat melakukan ini.

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

Mengevaluasi fungsi yang didefinisikan secara terpisah.

Dengan serangkaian kondisi dan fungsi terkait, evaluasi setiap fungsi pada input data di mana kondisinya benar.

Contoh diberikan pada SO di sini . Untuk kelengkapan, berikut adalah contohnya:

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

Metode yang diusulkan oleh Vito MR Muggeo [1] relatif sederhana dan efisien. Ini berfungsi untuk sejumlah segmen tertentu, dan untuk fungsi kontinu. Posisi breakpoint diperkirakan secara iteratif dengan melakukan, untuk setiap iterasi, regresi linier tersegmentasi yang memungkinkan lompatan pada breakpoint. Dari nilai-nilai lompatan, posisi breakpoint berikutnya disimpulkan, sampai tidak ada lagi diskontinuitas (melompat).

"prosesnya diulangi sampai kemungkinan konvergensi, yang secara umum tidak dijamin"

Secara khusus, konvergensi atau hasilnya mungkin tergantung pada estimasi pertama breakpoints.

Ini adalah metode yang digunakan dalam paket R Segmented .

Berikut ini adalah implementasi dalam 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

Contoh:

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

[1]: Muggeo, VM (2003). Memperkirakan model regresi dengan breakpoint yang tidak diketahui. Statistik dalam kedokteran, 22 (19), 3055-3071.

Saya sudah mencari hal yang sama, dan sayangnya sepertinya tidak ada saat ini. Beberapa saran untuk bagaimana melanjutkan dapat ditemukan dalam pertanyaan sebelumnya .

Atau Anda bisa melihat ke beberapa pustaka R misalnya tersegmentasi, SiZer, strucchange, dan jika ada sesuatu yang berfungsi untuk Anda, cobalah menyematkan kode R dalam python dengan rpy2 .

Mengedit untuk menambahkan tautan ke py-earth , "Sebuah implementasi Python dari Multivariate Adaptive Regression Splines Jerome Friedman".

Ada posting blog dengan implementasi regresi piecewise yang rekursif. Solusi itu cocok dengan regresi terputus-putus.

Jika Anda tidak puas dengan model terputus-putus dan ingin pengaturan terus menerus, saya akan mengusulkan untuk mencari kurva Anda dalam basis kkurva berbentuk-L, menggunakan Lasso untuk sparsity:

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

Kode ini akan mengembalikan vektor perkiraan koefisien kepada Anda:

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

Karena pendekatan Lasso, jarang: model menemukan tepat satu breakpoint di antara 10 kemungkinan. Angka 0,57 dan 0,825 sesuai dengan 0,5 dan 1,25 dalam DGP yang sebenarnya. Meskipun tidak terlalu dekat, kurva yang dipasang adalah:

Pendekatan ini tidak memungkinkan Anda memperkirakan breakpoint dengan tepat. Tetapi jika dataset Anda cukup besar, Anda bisa bermain dengan yang berbeda k(mungkin tune dengan cross-validation) dan perkirakan breakpoint dengan cukup tepat.