Last active
July 3, 2024 23:21
-
-
Save ruoyu0088/70effade57483355bbd18b31dc370f2a to your computer and use it in GitHub Desktop.
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
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?