10
$\begingroup$

Looking for a pleasing way to smooth a particular dataset I found something I like but it is rather slow to compute: I take the mean of a series of MedianFilter results computed over a range of window sizes.

SeedRandom[0];
x = RandomReal[1, 5000];

smooth = Mean @ Table[MedianFilter[x, n], {n, 500, 900, 5}];  (* slow *)

ListLinePlot[smooth]

enter image description here

(The closer the spacing the better the result to my eye, i.e. keeping 5 a small value.)

Is there a faster way to get the same smoothing, or something very similar?

Note:

Separate from this question is the issue of whether or not this smoothing actually represents the data in a useful or accurate way. I would rather not focus on that but instead accept this result being desirable as an assumption.

$\endgroup$
  • 1
    $\begingroup$ The sense I get is that you want a smoothing without any respect to an underlying model that generates such data. If so and speed is an issue why not use MeanFilter rather than MedianFilter? $\endgroup$ – JimB Dec 27 '16 at 6:43
  • 1
    $\begingroup$ I tossed out this question before going out of town wondering what I'd find when I came back. I am impressed with the attention it has been given and I thank you all for your time. I know the question is not well formed with a crummy example dataset and little explanation and I am sorry for that; had I known the effort people were going to give this I would have tried harder. I am still tired myself but I shall look at all these answers in more depth in the next few days. I'll also try to improve the question a bit without "moving the goalposts" if I am able. Thanks again. $\endgroup$ – Mr.Wizard Dec 31 '16 at 3:08
  • $\begingroup$ Documentation of NetTrain(ref/NetTrain-Applications-Manifold Learning) maybe can help though I cannot reproduce it. $\endgroup$ – yode Feb 15 '17 at 19:05
10
$\begingroup$

All smooths in this post

Here are all smooths in this post:

enter image description here

The line "smooth" is the one derived in the question. The lines for "smooth2" and "smooth3" were the motivation to derive the Quantile regression solution, "smoothQR". The speed-ups are around 7 and 4.5 times respectively. (See the related discussion in Rahul's answer.)

Simple modifications of question's code

Switching Mean and Median makes the computations ~ 7 times faster on my laptop. (The rationale for doing that is that finding the median is slower than finding the mean.)

SeedRandom[0];
x = RandomReal[1, 5000];

AbsoluteTiming[
 smooth = Mean@Table[MedianFilter[x, n], {n, 500, 900, 5}];
](* original code in the question *)

(* {8.44073, Null} *)

ListLinePlot[smooth]

enter image description here

Mean and Median switched around:

AbsoluteTiming[
 smooth2 = Median@Table[MeanFilter[x, n], {n, 500, 900, 5}];
]

(* {1.31114, Null} *)

ListLinePlot[smooth2]

enter image description here

Or using just Mean:

AbsoluteTiming[
 smooth3 = Mean@Table[MeanFilter[x, n], {n, 500, 900, 5}];
]

(* {1.41859, Null} *)

ListLinePlot[smooth3]

enter image description here

Using Quantile regression

Here is a more complicated way to speed up things together with some error measurements.

  1. Select parameters.

  2. Apply MedianFilter as in the question but with larger step between the min and max window sizes.

  3. Find the mean or median of the obtained array. (We get a vector as a result.)

  4. Apply quantile regression with B-spline basis over a sample of the points from step 3.

  5. Using the found regression quantile function extract points.

  6. Plot and measure errors.

The process of steps 2-5 is around 3 to 5 times faster on my laptop than the original code in the question and with a reasonable max relative error, less than 0.5%.

{a, b, step} = {500, 900, 5};
{factor, sampleSize, nKnots} = {16, 800, 25};

AbsoluteTiming[
 smooth = Mean@Table[MedianFilter[x, n], {n, a, b, s}];
]

(* {8.62565, Null} *)

Import["https://raw.githubusercontent.com/antononcube/MathematicaForPrediction/master/QuantileRegression.m"]

AbsoluteTiming[
 smoothArr = Table[MedianFilter[x, n], {n, a, b, factor*step}]; 
 ts = N@Transpose[{Rescale@Range[Length[Mean[smoothArr]]], 
     Mean[smoothArr]}]; 
 qfuncs = QuantileRegression[
   Join[ts[[1 ;; 10]], 
    RandomSample[ts[[21 ;; -11]], UpTo[sampleSize]], ts[[-10 ;; -1]]],
    nKnots, {0.5}, 
   Method -> {LinearProgramming, Method -> "InteriorPoint", 
     Tolerance -> 1.}, InterpolationOrder -> 3];
 smoothQR = 
  Table[qfuncs[[1]][y], {y, ts[[1, 1]], ts[[-1, 1]], (
    ts[[-1, 1]] - ts[[1, 1]])/(Length[x] - 1)}];
 ]

(* {1.72669, Null} *)

grOpts = {PlotRange -> All, PlotTheme -> "Detailed", 
  ImageSize -> 450}; ListLinePlot[{smooth, smoothQR}, 
 PlotLegends -> SwatchLegend[{"smooth", "smoothQR"}], grOpts]
ListPlot[Abs[smooth - smoothQR]/Abs[smooth], 
 PlotLegends -> 
  RecordsSummary[
   Abs[smooth - smoothQR]/Abs[smooth], {"Relative error"}], 
 PlotLabel -> 
  Row[{"Relative error, ", 
    HoldForm[Abs[smooth - smoothQR]/Abs[smooth]]}], grOpts]

enter image description here

$\endgroup$
7
$\begingroup$

One way to look at your method is that for each sample, the samples within distance 500 of it get considered in all of the medians you are taking the mean of, the samples at distance 700 get considered in half the medians, the samples at distance 900 get considered in only one, and so on.

What if we take just one median, but with each sample weighted according to the number of times it would have been considered in your original approach?

wmin = 500; 
wmax = 900; 
newSmooth = 
 Table[With[{range = Range[Max[i - wmax, 1], Min[i + wmax, Length@x]]}, 
   Median@WeightedData[x[[range]], 
     Min[#, wmax - wmin + 1] & /@ (wmax + 1 - Abs[range - i])]], 
   {i, Length@x}];

ListLinePlot[{smooth, newSmooth}, PlotStyle -> {Gray, Red}]

enter image description here

Pretty close, and about 10 times faster on my machine; switching to ParallelTable adds another 3x speedup. Though the result does have those telltale MedianFilter-esque discontinuities that you might not like.

(Note: The code above is designed so that if you choose wmin = wmax the result should be identical to MedianFilter[x, wmax], while if you set wmin = 0 you should get an exactly triangular weight distribution.)


Since Mr.Wizard has been cagey about their actual data and motivation (They should really know better by now. Maybe a moderator could knock some sense into them.), here is an example dataset that illustrates why one might prefer a MedianFilter-based approach.

x = ImageData[ExampleData[{"TestImage", "Flower"}]][[1500, All, 1]];

Here's Mr.Wizard's smoothing method, though I changed the smoothing range to 100-500 instead of 500-900:

Show[ListPlot[x, PlotStyle -> Lighter@Gray], 
 ListLinePlot[smooth, PlotStyle -> Black]]

enter image description here

As you can see, MedianFilter ignores narrow spikes of outliers without blurring across significant discontinuities. Here are the rest of the smoothing methods for comparison:

Show[ListPlot[x, PlotStyle -> LightGray], 
 ListLinePlot[{smooth, smooth2, smooth3, newSmooth}, 
  PlotStyle -> {Black, Darker@Cyan, Lighter@Blue, Lighter@Red}, 
  PlotLegends -> {"Mr.Wizard", "Anton2", "Anton3", "Rahul"}]]

enter image description here

$\endgroup$
  • $\begingroup$ "Maybe a moderator could knock some sense into them." - :D Sadly I cannot upvote again. $\endgroup$ – J. M. is in limbo Dec 28 '16 at 20:45
  • $\begingroup$ +1 for the comment about motivation (although +1 in any event). Data manipulation techniques without (stated) purpose are close to useless. That goes for choosing between means and medians. But note that spikes can be good data, too! Not all spikes are outliers! $\endgroup$ – JimB Dec 28 '16 at 20:53
  • $\begingroup$ @Jim: I agree! Though by definition all spikes are outliers; don't you mean that not all outliers are bad data? :) $\endgroup$ – Rahul Dec 28 '16 at 21:06
  • 1
    $\begingroup$ One first needs to define what an "inlier" before determining what an outlier is. An outlier requires a purpose that you've noted is missing. We have different definitions for spikes. In chromatograms of cuticular hydrocarbons a collection of spikes form a unique signature to identify termite species (as one of many examples) and smoothing just wipes out the signature. $\endgroup$ – JimB Dec 28 '16 at 21:15
  • $\begingroup$ @Jim: The fact that smoothing wipes out the signature need not be a bad thing, either. One can then subtract the smoothed signal from the original to get the spikes cleanly isolated from any slowly-varying background, e.g. the Matlab documentation on correcting nonuniform illumination. $\endgroup$ – Rahul Dec 28 '16 at 21:34
5
$\begingroup$

Only considering speeding it up you may use ParallelMap to spread the filtering over multiple cores.

LaunchKernels[];
smooth = ParallelMap[MedianFilter[x, #] &, Range[500, 900, 5]];
CloseKernels[];

smooth = Mean@smooth;

This executes in 13 seconds on my laptop whereas the non-parallel code takes 32 seconds.

Hope this helps.

$\endgroup$
5
$\begingroup$

Rather than an answer, this is an extended comment. All of the answers so far speed things up considerably. However, the size of the observed gyrations of the smoothed data needs to be considered relative to the scatter of points:

Show[ListPlot[x],
 ListLinePlot[smooth, PlotStyle -> {Red, Thickness[0.02]}]]

Point scatter and smoothed curve

The process that generated the example dataset (which I assume is much different than your real dataset of interest) should result in a very straight horizontal line with the vertical axis value of 0.5 with there being random but serially correlated deviations from that value.

The point is that some sort of regression approach (loess, quantile regression, generalized additive models, etc.) might be considered prior to smoothing to test that any of the observed "features" might be real (rather than expected random deviations).

$\endgroup$
3
$\begingroup$

Edit: after working this up I realized I've mimicked MeanFilter, not MedianFilter

SeedRandom[0];
x = RandomReal[1, 5000];
wt = Join[#, ConstantArray[81, 1000], Reverse@#] &@
   Flatten[Table[ ConstantArray[ j, 5] , {j, 80}]];
filter[x_, i_] := Module[{wl = Round[Length@wt/2]}, Which[ 
   i < wl , (x[[;; i + wl]].#/Total@#) &@wt[[-(wl + i) ;;]],
   i > Length@x - wl , (x[[i - wl ;;]].#/Total@#) &@
    wt[[;; Length@x - wl - i]],
   True, x[[i - Length@wt/2 + 1 ;; i + wl]].wt/(Total@wt)]]

ListLinePlot[filter[x, #] & /@ Range[Length@x]]

enter image description here

$\endgroup$
3
$\begingroup$

This is fast, and gives a similar smoothing, though probably not as similar as you would like. It relies on the data in x being uniformly distributed though I think it could be adapted for other cases.

Consider a run of $m$ elements from $x$. If there are $m/2$ points less than 0.5 and $m/2$ greater than 0.5 then the median will be somewhere very close to 0.5. Changing one of the <0.5 points to become >0.5 will shift the median by one place, on average an increase of $1/m$. So we can approximate the median by looking at the number of points either side of 0.5. Each point <0.5 shifts the median by $-1/(2m)$ and each point >0.5 shifts the median by $+1/(2m)$.

The result can be computed quickly using ListCorrelate. It runs in a few hundredths of a second.

y = Sign[x - 0.5];

ca[n_] := ConstantArray[0.5/(2 n + 1), 2 n + 1]

smooth2 = 
  Mean @ Table[0.5 + ListCorrelate[ca[n], y, {n + 1, n + 1}, 0], {n, 500, 900, 5}];

ListLinePlot[{smooth, smooth2}]

enter image description here

$\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.