18
$\begingroup$

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.

enter image description here

|improve this question|||||
$\endgroup$
  • $\begingroup$ I don't think this question can be answered reasonably unless you tell us how accurately you want to know the locations of the break-points, what your guesstimate is for the shortest length of a linear segment and how many samples there are in a typical transition region. If the horizontal axis labels in your figure are sample numbers, then, with two transitions in the span from $x[-5]$ to $x[0]$, the task is more difficult than if the straight-line segments were of longer duration (in samples). $\endgroup$ – Dilip Sarwate Jan 25 '12 at 16:54
  • $\begingroup$ @DilipSarwate I updated the Question with the requirements(btw the xaxis is the magnetic field in tesla) $\endgroup$ – P3trus Jan 25 '12 at 17:59
  • $\begingroup$ You can try this toolbox if you are working with MATLAB curve fitting toolbox $\endgroup$ – Rhei Jan 18 '15 at 17:21
12
$\begingroup$

I tried two approaches, naively (using only 3 segments). Surely there would be fancier methods out there.

    RANSAC, supposed to be a robust fitting mechanism. It's easy to stop the algorithm after a number of segments. However it may be difficult to enforce continuity between segments--as seems required in your application-- at least with a simple implementation. As a proof of concept, I created an image from the data points so that I could use the RANSAC engine available in $ImageLines$, the line detection function of Mathematica.

enter image description here

    Fit a piecewise linear model using a general purpose minimizer. It's easy to enforce segments continuity. Interestingly, testing for residuals and other properties may provide enough information to determine automatically the number of segments--I've not tried it though. That's how it looks in Mathematica:

enter image description here

|improve this answer|||||
$\endgroup$
  • $\begingroup$ Looks like a great answer. Thanks for contributing. $\endgroup$ – Jason R Jan 26 '12 at 13:51
7
$\begingroup$

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.

|improve this answer|||||
$\endgroup$
  • $\begingroup$ Totally stole my answer! = ) $\endgroup$ – Phonon Jan 26 '12 at 19:21
  • $\begingroup$ Interresting Idea but sadly due to noise on the signal I don't get good results. $\endgroup$ – P3trus Jan 28 '12 at 6:11
  • 1
    $\begingroup$ That expression whose magnitute is being compared to epsilon is actually an approximation to the second derivative of the data. There are other ways to compute this using more than three point which don't respond to noise as much. Look up Savitzky-Golay. $\endgroup$ – DarenW Nov 11 '12 at 2:32
4
$\begingroup$

(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.

Added March 2014: Dynamic programming is a general method for problems with nested subproblems like this:

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.

enter image description here

|improve this answer|||||
$\endgroup$
  • $\begingroup$ The problem, at least with scipy is the positioning of the knots. scipy uses equaly spaced knots. $\endgroup$ – P3trus Sep 17 '13 at 8:49
  • $\begingroup$ @P3trus, yes for a start, but then they can move -- see the plot. Anyway it targets total error, not knots. $\endgroup$ – denis Sep 18 '13 at 12:29
  • $\begingroup$ @P3trus Have you tried using multivariate regression splines method that automatically selects the breakpoints iteratively? cs.rtu.lv/jekabsons/regression.html $\endgroup$ – Atul Ingle Oct 13 '13 at 20:21
  • $\begingroup$ @Atul Ingle, afaik breakpoint / knot selection is the same problem, from whatever spline fitter. If you know of different algorithms for that from R / regression people, could you post a link please ? $\endgroup$ – denis Oct 25 '13 at 10:11
  • $\begingroup$ Are looking for packages in R/Matlab that do adaptive regression splines? Here: cran.r-project.org/web/packages/earth/index.html cran.r-project.org/web/packages/mda/index.html and also ARESLab in Matlab that I already posted the link for. $\endgroup$ – Atul Ingle Oct 30 '13 at 23:28
0
$\begingroup$

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.

|improve this answer|||||
$\endgroup$
  • $\begingroup$ derivative might be offul noisy. i don't think i would recommend that. $\endgroup$ – robert bristow-johnson Mar 10 '14 at 3:47
0
$\begingroup$

Using an l1 trend filter is another idea:

Paper

Online Example

|improve this answer|||||
$\endgroup$
  • 1
    $\begingroup$ Your answer is a little bit too short to be constructive! Please consider making an effort to expand it in a pedagogical way. $\endgroup$ – sansuiso Apr 30 '14 at 7:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.