segments_fit.ipynb

segments_fit2.ipynb

Thanks for posting this!

Here are a few suggested edits to pick the number of segments automatically by optimizing AIC and/ or BIC. The implementation is similar to the heuristic strategy presented in this paper: https://discovery.ucl.ac.uk/id/eprint/10070516/1/AIC_BIC_Paper.pdf

def segments_fit(X, Y, maxcount):
    xmin = X.min()
    xmax = X.max()
    
    n = len(X)
    
    AIC_ = float('inf')
    BIC_ = float('inf')
    r_   = None
    
    for count in range(1, maxcount+1):
        
        seg = np.full(count - 1, (xmax - xmin) / count)

        px_init = np.r_[np.r_[xmin, seg].cumsum(), xmax]
        py_init = np.array([Y[np.abs(X - x) < (xmax - xmin) * 0.1].mean() for x in px_init])

        def func(p):
            seg = p[:count - 1]
            py = p[count - 1:]
            px = np.r_[np.r_[xmin, seg].cumsum(), xmax]
            return px, py

        def err(p): # This is RSS / n
            px, py = func(p)
            Y2 = np.interp(X, px, py)
            return np.mean((Y - Y2)**2)

        r = optimize.minimize(err, x0=np.r_[seg, py_init], method='Nelder-Mead')
    
        # Compute AIC/ BIC. 
        AIC = n * np.log10(err(r.x)) + 4 * count
        BIC = n * np.log10(err(r.x)) + 2 * count * np.log(n)
        
        if((BIC < BIC_) & (AIC < AIC_)): # Continue adding complexity.
            r_ = r
            AIC_ = AIC
            BIC_ = BIC
        else: # Stop.
            count = count - 1
            break
        
    return func(r_.x) ## Return the last (n-1)

Hi dankoc! Thanks for sharing the code which really helps me. One little clarification, shouldn't it be natural log in the AIC/BIC formula?

Yes, it should use the natural log.

…

On Mon, Feb 20, 2023 at 4:58 PM Qinkong ***@***.***> wrote: ***@***.**** commented on this gist. ------------------------------ Thanks for posting this! Here are a few suggested edits to pick the number of segments automatically by optimizing AIC and/ or BIC. The implementation is similar to the heuristic strategy presented in this paper: https://discovery.ucl.ac.uk/id/eprint/10070516/1/AIC_BIC_Paper.pdf def segments_fit(X, Y, maxcount): xmin = X.min() xmax = X.max() n = len(X) AIC_ = float('inf') BIC_ = float('inf') r_ = None for count in range(1, maxcount+1): seg = np.full(count - 1, (xmax - xmin) / count) px_init = np.r_[np.r_[xmin, seg].cumsum(), xmax] py_init = np.array([Y[np.abs(X - x) < (xmax - xmin) * 0.1].mean() for x in px_init]) def func(p): seg = p[:count - 1] py = p[count - 1:] px = np.r_[np.r_[xmin, seg].cumsum(), xmax] return px, py def err(p): # This is RSS / n px, py = func(p) Y2 = np.interp(X, px, py) return np.mean((Y - Y2)**2) r = optimize.minimize(err, x0=np.r_[seg, py_init], method='Nelder-Mead') # Compute AIC/ BIC. AIC = n * np.log10(err(r.x)) + 4 * count BIC = n * np.log10(err(r.x)) + 2 * count * np.log(n) if((BIC < BIC_) & (AIC < AIC_)): # Continue adding complexity. r_ = r AIC_ = AIC BIC_ = BIC else: # Stop. count = count - 1 break return func(r_.x) ## Return the last (n-1) Thanks Thanks for posting this! Here are a few suggested edits to pick the number of segments automatically by optimizing AIC and/ or BIC. The implementation is similar to the heuristic strategy presented in this paper: https://discovery.ucl.ac.uk/id/eprint/10070516/1/AIC_BIC_Paper.pdf def segments_fit(X, Y, maxcount): xmin = X.min() xmax = X.max() n = len(X) AIC_ = float('inf') BIC_ = float('inf') r_ = None for count in range(1, maxcount+1): seg = np.full(count - 1, (xmax - xmin) / count) px_init = np.r_[np.r_[xmin, seg].cumsum(), xmax] py_init = np.array([Y[np.abs(X - x) < (xmax - xmin) * 0.1].mean() for x in px_init]) def func(p): seg = p[:count - 1] py = p[count - 1:] px = np.r_[np.r_[xmin, seg].cumsum(), xmax] return px, py def err(p): # This is RSS / n px, py = func(p) Y2 = np.interp(X, px, py) return np.mean((Y - Y2)**2) r = optimize.minimize(err, x0=np.r_[seg, py_init], method='Nelder-Mead') # Compute AIC/ BIC. AIC = n * np.log10(err(r.x)) + 4 * count BIC = n * np.log10(err(r.x)) + 2 * count * np.log(n) if((BIC < BIC_) & (AIC < AIC_)): # Continue adding complexity. r_ = r AIC_ = AIC BIC_ = BIC else: # Stop. count = count - 1 break return func(r_.x) ## Return the last (n-1) Hi dankoc! Thanks for sharing the code which really helps me. One little clarification, shouldn't it be natural log in the AIC/BIC formula? — Reply to this email directly, view it on GitHub <https://gist.github.com/70effade57483355bbd18b31dc370f2a#gistcomment-4477503> or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAYUH7OJLOUYQZRZ5DFMF6LWYPSH7BFKMF2HI4TJMJ2XIZLTSKBKK5TBNR2WLJDHNFZXJJDOMFWWLK3UNBZGKYLEL52HS4DFQKSXMYLMOVS2I5DSOVS2I3TBNVS3W5DIOJSWCZC7OBQXE5DJMNUXAYLOORPWCY3UNF3GS5DZVRZXKYTKMVRXIX3UPFYGLK2HNFZXIQ3PNVWWK3TUUZ2G64DJMNZZDAVEOR4XAZNEM5UXG5FFOZQWY5LFVA2DONZZGY4TGOFHORZGSZ3HMVZKMY3SMVQXIZI> . You are receiving this email because you commented on the thread. Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub> .

-- Associate Professor Baker Institute for Animal Health College of Veterinary Medicine Cornell University Website: http://www.dankolab.org E-mail: ***@***.***

ruoyu0088 Thanks for the original version - it was exactly what I needed for a piecewise linear approximation of a CDF I wanted.

First off, thanks a bunch for this, @ruoyu0088 ! One of the fastest solvers I've found outthere. One small hiccup though, as you can see, the segment 2 is going back in time, which is a big no-no in the timeseries world and I could not figure out how to fix that

@josesydor

Please try the following code:

from scipy import optimize

def segments_fit(X, Y, count):
    xmin = X.min()
    xmax = X.max()
    ymin = Y.min()
    ymax = Y.max()    

    seg = np.full(count - 1, (xmax - xmin) / count)

    px_init = np.r_[np.r_[xmin, seg].cumsum(), xmax]
    py_init = np.array([Y[np.abs(X - x) < (xmax - xmin) * 0.01].mean() for x in px_init])

    def func(p):
        seg = p[:count - 1]
        py = p[count - 1:]
        px = np.r_[np.r_[xmin, seg].cumsum(), xmax]
        return px, py

    def err(p):
        px, py = func(p)
        Y2 = np.interp(X, px, py)
        return np.mean((Y - Y2)**2)

    r = optimize.minimize(
        err, x0=np.r_[seg, py_init], method='SLSQP', 
        bounds=[(0, None)] * len(seg) + [(ymin, ymax)] * len(py_init), 
        constraints=[
            optimize.LinearConstraint([1] * len(seg) + [0] * len(py_init), 0, xmax - xmin)
        ])
    
    return func(r.x)

@ruoyu0088 Thanks, but unfortunatelly it didnt work:

These are my X, Y:
3;6.37
4;6.38
5;6.38
11;6.42
12;6.37
14;6.33
15;6.34
16;6.35
17;6.35
19;6.37
20;6.4
21;6.41
22;6.39
23;6.37
24;6.39
25;6.4
26;6.35
27;6.32
28;6.31
30;6.31
31;6.3
32;6.3
33;6.29
34;6.29
35;6.27
36;6.28
37;6.28
38;6.29
40;6.34
41;6.33
42;6.33
43;6.35
45;6.36
46;6.35
47;6.36
49;6.33
50;6.35
51;6.35
52;6.35
54;6.35
55;6.34
56;6.34
57;6.35
58;6.35
59;6.34
60;6.33
61;6.32
62;6.33
63;6.32
64;6.31
65;6.31
66;6.31
67;6.31
69;6.32
70;6.3
72;6.29
73;6.29
74;6.29
75;6.3
76;6.28
77;6.28
78;6.29
79;6.29
80;6.3
81;6.31
82;6.32
84;6.32
86;6.33
87;6.33
88;6.35
90;6.33
92;6.34
93;6.33
94;6.32
96;6.33
97;6.33
98;6.34
99;6.33
101;6.31
102;6.32
103;6.32
105;6.32
106;6.33
108;6.33
109;6.31
110;6.31
111;6.34
112;6.35
114;6.35
116;6.34
117;6.35
118;6.35
119;6.34
120;6.34
121;6.33
123;6.32
124;6.32
125;6.32
126;6.33
128;6.32
129;6.33
130;6.33
131;6.32
132;6.33
133;6.33
134;6.32
135;6.32
136;6.32
137;6.32
138;6.32
139;6.32
141;6.32
143;6.33
144;6.33
145;6.32
146;6.32
147;6.32
148;6.3
149;6.28
150;6.29
151;6.28
152;6.28
153;6.27
154;6.26
155;6.26
156;6.28
157;6.25
158;6.25
159;6.24
160;6.24
161;6.22
162;6.21
163;6.2
164;6.18
165;6.18
166;6.18
168;6.17
169;6.17
170;6.16
171;6.15
172;6.15
173;6.18
174;6.16
175;6.16
176;6.16
177;6.17
178;6.15
179;6.15
180;6.16
181;6.16
182;6.16
183;6.16
184;6.16
185;6.16
186;6.17
187;6.16
188;6.16
189;6.16
190;6.16
191;6.16
192;6.16
193;6.17
194;6.16
195;6.15
196;6.16
197;6.15
198;6.16
199;6.14
200;6.14
201;6.13
202;6.13
203;6.13
204;6.13
205;6.13
206;6.13
207;6.14
208;6.14
209;6.14
210;6.14
211;6.14
212;6.14
213;6.15
214;6.15
215;6.14
216;6.12
217;6.12
218;6.12
219;6.12
220;6.11
221;6.12
222;6.11
223;6.11
224;6.09
225;6.11
226;6.1
227;6.09
228;6.1
229;6.1
230;6.1
231;6.1
232;6.09
233;6.09
235;6.1
237;6.08
238;6.08
239;6.08
240;6.09
241;6.08
242;6.1
243;6.11
244;6.1
245;6.12
246;6.13
247;6.14
248;6.15
250;6.13
251;6.12
252;6.11
253;6.12
254;6.14
255;6.14
256;6.15
257;6.13
258;6.12
259;6.1
260;6.12
261;6.12
262;6.12
263;6.11
264;6.12
265;6.13
266;6.13
267;6.15
268;6.15
269;6.14
270;6.16
271;6.12
272;6.12
273;6.11
274;6.11
275;6.11
276;6.1
277;6.09
278;6.08
279;6.06
280;6.07
281;6.07
282;6.08
283;6.08
284;6.06
285;6.08
286;6.08
287;6.07
288;6.08
289;6.07
290;6.08
291;6.08
292;6.08
293;6.09
294;6.08
295;6.08
296;6.06
297;6.05
298;6.05
299;6.06
300;6.05
301;6.06
302;6.05
303;6.04
304;6.05
305;6.05
306;6.06
307;6.05
308;6.05
309;6.04
310;6.05
311;6.04
312;6.05
313;6.05
314;6.07
315;6.06
316;6.06
317;6.06
318;6.05
319;6.05
320;6.05

This is the function original version. Again, segment 2 is going backwards:

@josesydor I uploaded another notebook to process your data, it looks ok.