Velocity estimation by piecewise linear regression #569
Conversation
Start with simple OLS linear regression if no breakpoints requested Set up model dataclass
| diff_term = x - breakpoints[:, np.newaxis] | ||
| heavi = np.vstack([np.heaviside(d, -1) for d in diff_term]) | ||
| u = diff_term * heavi | ||
| v = -heavi | ||
| design_matrix = np.vstack((np.ones(x.size), x, u, v)).T | ||
|
|
||
| xtx = np.linalg.pinv(np.matmul(design_matrix.T, design_matrix)) | ||
| coeffs = np.dot( | ||
| xtx, | ||
| np.dot(design_matrix.T, y), | ||
| ) | ||
|
|
||
| fit = np.dot(design_matrix, coeffs) |
There was a problem hiding this comment.
I wonder if moving the linear fit part and the parameter remapping inside PiecewiseModel would make sense. Since it's a linear fit, the fit is fully determined from just the breakpoints and the data (which could be its constructor arguments).
The benefit would be that you could add a simple .plot() function to that model to see it (since you'd have all the pieces there to plot it already).
Then you have a self-contained PiecewiseModel model with everything needed to fit (which can just happen on construction making it immutable) and assess its quality for a given set of breakpoints in a single dataclass.
bic, prediction, rss, breaks_in_range and external parameters could be properties.
The only thing that seems like it wouldn't fit quite as nicely is the exitflag; but that feels a bit like a property of the overarching algorithm rather than the PiecewiseLinearModel itself.
| from a uniform distribution. Ignored if `n_breakpoints == 0`. | ||
| """ | ||
| n_samples = len(x) | ||
| n_coeffs = 2 + 2 * n_breakpoints |
There was a problem hiding this comment.
I do wonder if BIC treating the breakpoint as just another parameter is sufficient penalty. In a sense, it is a very different type of parameter. See also the answer on this for example: https://stats.stackexchange.com/questions/337852/aic-bic-for-a-segmented-regression-model
When you run this with many simulations, how consistently does the BIC selection procedure recover the correct model?
| v = -heavi | ||
| design_matrix = np.vstack((np.ones(x.size), x, u, v)).T | ||
|
|
||
| xtx = np.linalg.pinv(np.matmul(design_matrix.T, design_matrix)) |
There was a problem hiding this comment.
There's a shorthand for matmul.
| xtx = np.linalg.pinv(np.matmul(design_matrix.T, design_matrix)) | |
| xtx = np.linalg.pinv(design_matrix.T @ design_matrix) |
| slope_terms = np.hstack((alpha, beta)) | ||
| slope_block = cov[1 : n_terms + 1, 1 : n_terms + 1] | ||
| slopes = np.array([np.sum(slope_terms[:j]) for j in range(1, n_terms + 1)]) | ||
| slopes_std = np.array([np.sqrt(np.sum(slope_block[:j, :j])) for j in range(1, n_terms + 1)]) |
There was a problem hiding this comment.
How close are these estimates to say a bootstrap with different realizations?
| break | ||
|
|
||
| # update breakpoints and check in bounds | ||
| breakpoints = gamma / beta + breakpoints |
There was a problem hiding this comment.
Any reason why we're not just doing:
| breakpoints = gamma / beta + breakpoints | |
| breakpoints += gamma / beta |
work in progress