-
Notifications
You must be signed in to change notification settings - Fork 0
/
segment.py
51 lines (36 loc) · 1.48 KB
/
segment.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
# breakpoint analysis
# copied from https://datascience.stackexchange.com/questions/8457/python-library-for-segmented-regression-a-k-a-piecewise-regression
# accessed Jun 2019
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