Last active
July 3, 2024 23:21
-
-
Save ruoyu0088/70effade57483355bbd18b31dc370f2a to your computer and use it in GitHub Desktop.
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: ***@***.***
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Hello, thanks for the material, really useful! I have a question: is there any such code available to accomodate for 3 dimensions?