23
$\begingroup$

I'm looking at this sequence of values:

enter image description here

I'd like to detect the points where the center of the time-series shifts (around x=1000 and x=2000). Many of the transforms and smoothing methods I have tried destroy this information.

We can visually see that there are 3 "major" subsequences centered around y=400, 200, and then 400 again, but I'm not sure how programmatically locate these major components as continuous subsequences.

How can I automatically detect what the underlying linear step function is here?

To get going, you can grab the data like this:

data = Uncompress[FromCharacterCode[
  Flatten[ImageData[Import["http://i.stack.imgur.com/1D962.png"],"Byte"]]]]
|improve this question
$\endgroup$
  • $\begingroup$ I find this an interesting question. You are looking for cluster points in the series. Are there more than two? Do you know them in advance? (In this case it would be not so hard.) Maybe it's not the right place here to ask - hence no answers. Asking in a math forum? $\endgroup$ – Darwin1871 Jan 28 '16 at 0:55
  • $\begingroup$ Maybe you could do a windowed mean and standard deviation and make a hypothesis test whether the next (or the +k-th) data point lies within a certain confidence interval. Doing this with different window sizes. $\endgroup$ – Darwin1871 Jan 28 '16 at 1:12
  • $\begingroup$ Why not a Bai-Perron test for structural breaks? Perron has Matlab code on his website which is very easy to use. $\endgroup$ – Titus Feb 11 at 20:55
33
$\begingroup$

I had a go with HiddenMarkovProcess[], based on the assumption that the data is normally distributed around two different means (it looks like it!). This approach should be fine for cases where the number of "states" is small, e.g. 2 in this case. Otherwise you're looking at Infinite Hidden Markov Models, or see the bottom of this answer.

To remove some spurious detections, I first applied a median filter to smooth the data. You could also chop off the first and last 50 points (that drop to zero) to improve some of the estimates:

(* data is the provided tabulated list *)
ydata = MedianFilter[data[[All, 2]], 40];

hmm = EstimatedProcess[ydata, HiddenMarkovProcess[2, "Gaussian"]];
foundStates = 
  FindHiddenMarkovStates[ydata, hmm, "PosteriorDecoding"];

(* Extract the mean positions from the Markov model *)
hmmMeans = First@(# /. NormalDistribution -> List) & /@ Last@hmm;
(* {184.383, 391.369} *)

(* Now generate the piecewise data *)
meanFoundStates = foundStates /. Table[i -> hmmMeans[[i]], {i, 2}];

(* Now plot for comparison *)
ListLinePlot[{data[[All, 2]], meanFoundStates}]

enter image description here

Finally, you can detect the positions of the shift (e.g. at x=1000) by using:

FoldList[Plus, 1, Length@# & /@ Split[foundStates, #2 - #1 == 0 &]]
(* {1, 253, 907, 2044, 2946, 3143} *)

You can see that the "big" changes you refer to at x=1000 and x=2000 are actually picked up at 907 and 2044.

Another good alternative would be to make the most of RLink and use some of the packages for changepoint detection available in R. There are a few examples in this blog post, and I would recommend looking at:

  • bcp - Bayesian Change Point detection
  • ecp - Nonparametric Multiple Change Point Analysis
|improve this answer
$\endgroup$
  • 2
    $\begingroup$ That is cool... +1 $\endgroup$ – garej Jan 27 '16 at 10:42
  • $\begingroup$ Nice work, but I'm still wondering how to find the x points 1k and 2k where the boundaries are? $\endgroup$ – M.R. Jan 27 '16 at 16:08
  • $\begingroup$ @M.R. see my update $\endgroup$ – dr.blochwave Jan 27 '16 at 18:39
  • $\begingroup$ Is there a way to flatten out the discontinuities at {517, 539, 596, 632} and at {2218, 2237, 2368, 2392, 2708, 2735} so that the graph looks more like @dr.belisarius's graph? And in general there many be many more plateaus not just the two means, but n means, will your approach work for that too? $\endgroup$ – M.R. Jan 27 '16 at 21:29
  • $\begingroup$ @M.R. It will work if n is finite (otherwise you need to look into infinite Markov models). Have you looked into the R examples I posted above? They are very well-explained in the documentation and I've had great success with both. $\endgroup$ – dr.blochwave Jan 27 '16 at 21:32
18
$\begingroup$ 
ListPlot@{l1, msf = MeanShiftFilter[l1, IntegerPart[Length@l1/10], MedianDeviation@l1, 
                             MaxIterations -> 10]}

Mathematica graphics

And here are the detected means (assuming there are three):

fc = FindClusters[msf];
Mean /@ fc
( *{3.77282, 220.788, 387.444} *)
|improve this answer
$\endgroup$
  • $\begingroup$ This is a beautiful answer! How to I get the x values where it jumps? $\endgroup$ – M.R. Jan 27 '16 at 18:19
  • $\begingroup$ How and how did you know to divide by 10? $\endgroup$ – M.R. Jan 27 '16 at 18:20
  • $\begingroup$ It's perhaps worth pointing out that this doesn't actually give discrete values for each "state" - just a smoother version of the data. Using my method for detecting the jumps (FoldList[Plus, 1, Length@# & /@ Split[msf, #2 - #1 == 0 &]]) reveals this. $\endgroup$ – dr.blochwave Jan 27 '16 at 18:42
  • $\begingroup$ Although subsequent clustering around the means might work! $\endgroup$ – dr.blochwave Jan 27 '16 at 18:44
12
$\begingroup$

Another approach is to use compound median filtering which returns a blocky function. Then threshold the jumps between blocks. No assumptions about the number or size of blocks is made.

Function to plot the input series as discrete jumps.

BlockPlot[s_] := 
   Partition[
      Flatten[{s[[1]], 
         Table[{{s[[i, 1]], s[[i - 1, 2]]}, s[[i]]}, {i, 2, Length[s]}]}], 
   2]

Median filter until the signal does not change, then repeat for successively wider window radii r.

MedianFilterRoot[x_, r_] := FixedPoint[Round[MedianFilter[#, r]]&, x]

CompoundMedianFilter[x_?VectorQ, r_] := Fold[MedianFilterRoot[#1,#2]&,x,Range[r]]

CompoundMedianFilter[x_?MatrixQ, r_] :=
   Transpose[{x[[All,1]], Fold[MedianFilterRoot[#1, #2]&,x[[All,2]],Range[r]]}]

Find locations where the signal jumps more than the threshold t.

DifferenceThreshold[y_List, t_] :=
   Pick[Most[y], UnitStep[Abs[Differences[y[[All, 2]]]] - t], 1]

Subsample the imported data v down to w, for faster execution, and find the maximum. Both are global variables.

w = v[[Range[1, 3142, 10]]];
max = Max[w[[All, 2]]];

Adjust the maximum filter window radius r, and the jump threshold t.

Manipulate[
   Module[{y = CompoundMedianFilter[w, r]},
      ListLinePlot[BlockPlot[y], PlotStyle -> Directive[Thick, Blue],
         Prolog -> {Red, Point[w]},
         Epilog -> {Darker[Green],
            Map[Line[{{#[[1]], 0}, {#[[1]], max}}]&, DifferenceThreshold[y, t]]},
         Frame -> True, PlotRange -> {0, max}, ImageSize -> 600]],
   {{r, 10, "Max CMF Radius"}, 1, 30, 1, Appearance -> "Labeled"},
   {{t, 30, "Jump Threshold"}, 10, 200, 10, Appearance -> "Labeled"}
]

cmf

|improve this answer
$\endgroup$

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.