What is a robust way to fit piecewise linear but noisy data?
I'm measuring a signal, which consists of several almost linear segments. I'd like to atomatically fit several lines to the data to detect the transitions.
The dataset consists of a few thousand points, with 1-10 segments and I know the number of segments.
This is an example of what I'd like to do automatically.
I tried two approaches, naively (using only 3 segments). Surely there would be fancier methods out there.
I don't claim the following method is robust, but it might work for you. With thousands of points $x[n]$ and perhaps ten or so straight-line segments, proceed as follows.
Process the points $x[n]$ to create a bit array $y[n]$ as follows. $$y[n] = \begin{cases}1, &\text{if} ~ |(x[n+1]-x[n]) - (x[n]-x[n-1])| < \epsilon,\\ 0, &\text{otherwise.}\end{cases}$$ Here $\epsilon$ is a small number chosen to suit your notion of how close to a straight line you want points $x[n-1],x[n], x[n+1]$ to hew to. The criterion will be recognized by the cognoscenti as demanding that the straight line through $(n-1, x[n-1])$ and $(n,x[n])$ has nearly the same slope as the straight line through $(n,x[n])$ and $(n+1,x[n+1])$.
If $y[n]$ is an array of ten or so longish runs of $1$s separated by runs of $0$s with occasional stray $1$s here and there to mar the beauty, relax, you are on the right track. Else, if there are too few runs or too many runs of $1$s, repeat the previous step with a different $\epsilon$.
Use linear least-mean-square-error curve-fitting to fit straight lines to the points identified by $y[n]$ as belonging to the same straight-line segment. You now have ten straight lines fitting points, say, Line A fits points $x[3]$ through $x[88]$; line B fits points $x[94]$ through $x[120]$, Line C fits points $x[129]$ through $\cdots$, and so on. Extend A rightwards and B leftwards to find out where they intersect; extend B rightwards and C leftwards to find out where they intersect, etc. Congratulations, you now have a continuous and piecewise linear model for your data.
(Years later)
piecewise-linear functions are splines of degree 1, which most spline fitters can be told to do.
scipy.interpolate.UnivariateSpline
for example can be run with k=1
and a smoothing parameter s,
which you'll have to play with -- see
scipy-interpolation-with-univariate-splines .
In Matlab, see
how-to-choose-knots .
Added: finding optimal knots is not easy, because there can be many local optima.
Instead, you give UnivariateSpline a target s, sum of error^2,
and let it determine the number of knots.
After fitting, get_residual() will get the actual sum of error^2,
and get_knots() the knots.
A small change in s may change the knots a lot, especially in high noise -- ymmv.
The plot shows fits to a random piecewise-linear function + noise for various s.
For fitting piecewise constants, see Step detection. Can that be used for pw linear ? Don't know; starting off by differentiating noisy data will increase the noise, wrong.
Other testfunctions, and/or links to papers or code, would be welcome.
A couple of links:
piecewise-linear-regression-with-knots-as-parameters
$\qquad$ Linear splines are very sensitive to where the knots are placed
knot-selection-for-cubic-regression-splines
$\qquad$ This is a tricky problem and most people just select the knots by trial and error.
$\qquad$ One approach which is growing in popularity is to use penalized regression splines instead.
optimal k lines
= optimal k - 1 lines up to some x
+ cost of the last line x to the end
over x (all x in theory, nearby x in practice)
Dynamic programming is very clever, but can it beat brute force + heuristics for this task ?
See the excellent course notes by Erik Demaine under
MIT 6.006 Intro to algorithms
also google segmented linear regression
also John Henry syndrome.
Take the derivative and look for areas of nearly constant value. You would need to create the algorithm to search for those areas with ideally some level of +/- slope and that would give you the slope of the line for that section. You might want to perform some smoothing, like a sliding mean, prior to doing the sectional classification. The next step would be to get the y-intersection, which should be trivial at that point.
