segments_fit.ipynb

Hi! Sometimes, depending on the number of segments and the signals I pass, I get None in the return value. I believe it stems from the py_init initializing some indices to NaN. Why could it be? It leads me to a question: why do you multiply by 0.01?
Thank you!

@KateZi the multiplication factor of 0.01 ensures that only a small set of y-axis data (here 1% around the initial x-values) is selected to initialize the y-values corresponding to the initial x-values.

I had a data set which had a relatively larger difference between x-intervals than 1%, wherein I had to increase this parameter value to 0.02 (2%). That ensured that there are some y-values to initialize the optimization parameters.

Hello, thanks for the material, really useful! I have a question: is there any such code available to accomodate for 3 dimensions?

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: ***@***.***