6
$\begingroup$

I have the following data set:

  {{2, 25.4359}, {4, 25.1367}, {6, 24.8248}, {8, 24.8872},  
   {10, 24.7134}, {12, 24.3019}, {14, 23.2648}, {16, 22.7794}, 
   {18, 22.0979}, {20, 10.3315}, {22, 9.7597}, {24, 20.087}, 
   {26, 7.748}, {28, 17.3949}, {30, 5.952}, {32, 16.0672}, 
   {34, 15.084}, {36, 13.8762}, {38, 12.9525}, {40, 12.0011}, 
   {42, 11.1282}, {44, 10.1519}, {46, 9.56525}, {48, 8.80683}, 
   {50, 8.18156}, {52, 7.02375}, {54, 6.96017}, {56, 6.20118}, 
   {58, 5.57019}, {60, 5.82191}, {62, 5.61341}, {64, 5.40499}, 
   {66, 5.22668}, {68, 5.27839}, {70, 5.83658}}

enter image description here

Some of the data (red) in the range 20 < x 30 are too small (between 5 and 10) and have to be corrected to fit with the other data (blue).

How can this be achieved, without manually selecting the misplaced data?

I found that ClusterClassify can be used to distinguish the "wrong" data from the rest, but it could be accidental:

cl = ClusterClassify[data, 2, Method -> "Spectral"];

ListPlot[Pick[data, cl[data], #] & /@ {1, 2}, 
 FrameLabel -> {"x", "y"}, Frame -> True, PlotStyle -> {Blue, Red} ]

enter image description here

$\endgroup$
  • $\begingroup$ ...have to be corrected... If you intend to adjust the datapoints, you should have an idea of the error you wish to correct. From what I can see, the outliers are shifted downwards by a constant. Is that the case? Besides that, the use of ClusterClassify looks good, but it is likely to stumble on more difficult cases. $\endgroup$ – LLlAMnYP Feb 21 '18 at 17:13
  • $\begingroup$ @LLlAMnYP: One solution is probably after data separation to fit the blue data points and then to shift the red data points vertically upward so that they lie on the fitting curve. $\endgroup$ – mrz Feb 21 '18 at 17:39
  • 1
    $\begingroup$ In that case f = Interpolation[blue]; shift = (f /@ red[[;;,1]] - red[[;;,2]] // Mean); red[[;;,2]] += shift; data = Join[red,blue]//SortBy[First]. To move forward with an answer to your question a decent set of training data is necessary to optimize the parameters for ClusterClassify or a similar function that will detect outliers. $\endgroup$ – LLlAMnYP Feb 21 '18 at 17:59
  • $\begingroup$ For the specific dataset you could fit (as in regression) a sigmoid or some other relatively well behaved function and replace observations 'far' from the mean, with the mean; you should decide how far is 'far' eg one deviation. $\endgroup$ – yosimitsu kodanuri Feb 21 '18 at 17:59
  • 3
    $\begingroup$ Just to expand on the point made by @LLlAMnYP...Any adjustment method that doesn't acknowledge or account for the actual process that produces the points to be corrected is likely to fail. One needs to define both what an "outlier" is and an "inlier". Here is a recent but more complicated example (that is still unanswered): stats.stackexchange.com/questions/329170/…. $\endgroup$ – JimB Feb 21 '18 at 18:21
9
$\begingroup$

The procedure below finds the outliers, excludes them, fits the good data, then adds in the outliers corrected to fit the good data. The best fit appears to be a simple cubic function although I did try a sigmoid: y = c + (d - c)/(1 + Exp[a + b x]).

Outliers can be selected with variations on "stddev" -> 2 or "stddev" -> {2} to select outside or within band 2 respectively. The bands are 1, 2 & 3 standard deviations from the mean.

data = {{2, 25.4359}, {4, 25.1367}, {6, 24.8248}, {8, 24.8872},
   {10, 24.7134}, {12, 24.3019}, {14, 23.2648}, {16, 22.7794},
   {18, 22.0979}, {20, 10.3315}, {22, 9.7597}, {24, 20.087},
   {26, 7.748}, {28, 17.3949}, {30, 5.952}, {32, 16.0672},
   {34, 15.084}, {36, 13.8762}, {38, 12.9525}, {40, 12.0011},
   {42, 11.1282}, {44, 10.1519}, {46, 9.56525}, {48, 8.80683},
   {50, 8.18156}, {52, 7.02375}, {54, 6.96017}, {56, 6.20118},
   {58, 5.57019}, {60, 5.82191}, {62, 5.61341}, {64, 5.40499},
   {66, 5.22668}, {68, 5.27839}, {70, 5.83658}};

ListPlot[data]

enter image description here

outliers[data_, OptionsPattern["stddev" -> 2]] := Module[{},

   {sd1, sd2, sd3} = Map[2 (CDF[NormalDistribution[0, 1], #] - 0.5) &,
     {1, 2, 3}];

   Clear[x];

   lm = LinearModelFit[data, {1, x, x^2, x^3}, x];

   getseries[sd_] := Module[{cb},
     cb = lm["SinglePredictionBands", ConfidenceLevel -> sd];
     {lower, upper} = Transpose[cb /. x -> # & /@ First /@ data];
     p1 = Position[MapThread[#1 <= #2 <= #3 &,
        {lower, Last /@ data, upper}], True];
     Extract[data, p1]];

   If[ListQ[sd = OptionValue["stddev"]],
    o1 = Switch[sd,
      {1}, getseries[sd1],
      {2}, Complement[getseries[sd2], getseries[sd1]],
      {3}, Complement[getseries[sd3], getseries[sd2]],
      _, Print["\"stddev\" -> {2} used"];
      Complement[getseries[sd2], getseries[sd1]]],
    If[sd === None, o1 = data,
     s1 = getseries[Part[{sd1, sd2, sd3},
        Switch[sd, 1, 1, 2, 2, 3, 3,
         _, Print["\"stddev\" -> 2 used"]; 2]]];
     o1 = Complement[data, s1]]];

   {bands68[x_], bands95[x_], bands99[x_]} =
    Table[lm["SinglePredictionBands",
      ConfidenceLevel -> cl], {cl, {sd1, sd2, sd3}}];

   {{xmin, xmax}, {ymin, ymax}} = Through[{Min, Max}@#] & /@
     Transpose[data];

   scale[min_, max_, a_] := {(1 + a) min - a max, (1 + a) max + a min};

   {xmin1, xmax1} = scale[xmin, xmax, 0.05];
   {ymin1, ymax1} = scale[ymin, ymax, 0.2];
   {xmin2, xmax2} = scale[xmin, xmax, 0.1];

   Print[Show[ListPlot[data],
     Plot[{lm[x], bands68[x], bands95[x], bands99[x]},
      {x, xmin1, xmax1},
      Filling -> {2 -> {1}, 3 -> {2}, 4 -> {3}}],
     If[o1 != {}, ListPlot[o1,
       PlotStyle -> Directive[Blue, PointSize[Large]]], {}],
     AxesOrigin -> {xmin2, ymin1},
     PlotRange -> {{xmin2, xmax2}, {ymin1, ymax1}},
     ImageSize -> 400, Frame -> True]]];

outliers[data, "stddev" -> {3}]

enter image description here

data = Complement[data, o2 = o1];
outliers[data, "stddev" -> None]

enter image description here

o3 = Transpose[{First /@ o2,
    lm["BestFit"] /. # & /@ Thread[x -> First /@ o2]}];

Print[Show[Plot[{lm[x], bands99[x]},
   {x, xmin1, xmax1}, Filling -> {2 -> {1}}],
  ListPlot[o3,
   PlotStyle -> Directive[Red, PointSize[Large]]],
  ListPlot[data,
   PlotStyle -> Directive[Blue, PointSize[Large]]],
  AxesOrigin -> {xmin2, ymin1},
  PlotRange -> {{xmin2, xmax2}, {ymin1, ymax1}},
  ImageSize -> 400, Frame -> True]]

enter image description here

$\endgroup$
  • $\begingroup$ Thank you very much for your help. $\endgroup$ – mrz Feb 26 '18 at 6:50
4
$\begingroup$

An approach where each point is scanned and kept only if it is close enough from the previous one. The correction is done through an interpolating function on the good points.

threshold = 5;
i = 1; dataA = data;
l = Length@dataA;
While[i < l, 
 If[Abs[dataA[[i + 1, 2]] - dataA[[i, 2]]] < threshold
    , i++
    , dataA = Drop[dataA, {i + 1}]; l = Length@dataA
 ]
]
dataB = Complement[data, dataA];
if = Interpolation[dataA];
dataB2 = {#, if[#]} & /@ dataB[[All, 1]];
ListPlot[{dataA, dataB, dataB2}]

enter image description here

dataA is the list of good points, dataB the list of bad points, dataB2 the list of corrected points.

$\endgroup$
4
$\begingroup$

interpolate median filtered data:

medinterp = 
  Interpolation[
   Transpose@{data[[All, 1]], MedianFilter[data[[All, 2]], 6]}];

now drop the furthest outlier from the interpolation and redo the interpolation:

data = Drop[data, 
   First@Position[data, 
     Last@SortBy[data, Abs[medinterp[#[[1]]] - #[[2]]] &]]];
medinterp = 
  Interpolation[
   Transpose@{data[[All, 1]], MedianFilter[data[[All, 2]], 6]}];
Show[{
  Plot[medinterp[x], {x, 2, 70}],
  ListPlot[data]}]

manually eval the above 4 times..

enter image description here

fill in the points we dropped based on the interpolation:

 ListPlot[Join[
   data, {#[[1]], medinterp[#[[1]]]} & /@ Complement[data0, data]]]

enter image description here

note for this example you might do the same using a model fit (to Erf maybe)

$\endgroup$
2
$\begingroup$

One simple approach is to use the Median filter:

med = MedianFilter[data[[All, 2]], 4]
ListPlot[Transpose@{data[[All, 1]], med}]

enter image description here

$\endgroup$
  • $\begingroup$ Thank you. The problem is that I don't want to change the "blue" data points in my example (they represent correct measurements), but only "correct" the red ones so that they "fit" with the blue data. $\endgroup$ – mrz Feb 21 '18 at 19:30
  • $\begingroup$ @mrz -- It is a naive view of data which imagines that some data points are exactly correct and others are completely wrong. All measurements are inaccurate though maybe by different degrees, $\endgroup$ – bill s Feb 27 '18 at 1:00
2
$\begingroup$

Since there is no mention about the underlying data generating process, this answer is going to provide a simple way of distinguishing between outliers and regular data.

We'll assume that data (please evaluate the original expression before continuing) can be sensibly described by a sigmoid function like $f(x)=\frac{L}{1+e^{-k(x-x_0)}}$ which can be written also as $log(\frac{L}{f}-1)=kx_0-kx$. Note, how the last relation can be estimated using ordinary least squares only if $L$ is known.

Without any apriori knowledge about $L$ (besides the assertion in Wikipedia that "...L = the curve's maximum value" ) we will allow the user to set the appropriate value of $L$ according to her/his judgement.

Namely the code will allow the user to influence the value of $L$ by varying pL in (1 + pL ) ymax (the percentage by which $L$ exceeds the maximum value in data).

The 'knobs' in Manipulate control pL (f) and the radius within which points are considered non-outliers (r). Higher values for r allow more distant points from the estimated sigmoid to be considered normal.

Manipulate[

 display[pL, r],

 {{pL, 0.008, "f"}, 10^-9, 0.2, 10^-6},
 {{r, 2, "r"}, 1, 11, 0.5},

 Initialization :> (

   With[{xs = Part[data, All, 1], y = Part[data, All, -1], Nn = Length[data]},
    With[{ymax = Max[y]},

      (* ests will be used to construct the estimated data *)
      SetAttributes[ests, Listable];

      (* the sigmoid function that is assumed to describe the data *)
      s[x_] := L/(1 + Exp[-k (x - xo)]);

      (* Manipulate's first argument *)
      display[pL_, r_] := Module[{lmf, sol, dataest, time, junk, x, points, mask, n},
        (* transform data *)
        With[{dt = Transpose[{xs, Log[(1 + pL ) ymax/y - 1]}]},

          (* ols on transformed data *)
          lmf = LinearModelFit[dt, x, x];

          (* retrieve point estimates of parameters k and xo *)
          {sol} = Quiet[Solve[Thread[lmf["BestFitParameters"] == {k xo, -k}]]];
          (* augment solution with the value of L provided by the user *)
          AppendTo[sol, L -> (1 + pL) ymax];

          (* function used to produce the estimated data *)
          ests[x_] := Evaluate[s[x] /. sol];

          (* estimated data *)
          dataest = Transpose[{xs, ests[xs]}];

          (* indicates those elements in data that are further away from dataest than r ie the outliers *)
          mask = MapThread[EuclideanDistance[#1, #2] > r &, {dataest, data}];
          (* the actual points that correct for the outliers *)
          points = Pick[dataest, mask, True];

          n = Length[points];

          If[
             (* provided there are outliers *)
             n > 0,

             ListPlot[
               (* data excluding outliers, outliers, repaired data *)
               {Pick[data, Not /@ mask, True], Pick[data, mask, True], dataest},
               Joined -> {False, False, True},
               PlotStyle -> {
                 {PointSize[Medium]},
                 {PointSize[Medium], Opacity[0.4], Gray},
                 {Thin, Dashed, Opacity[0.6], Red}
                },
               PlotLegends -> 
                 Column[{Row[{"f=", N[pL, {2, 9}]}], 
                   Row[{"r=", N[r, {2, 9}]}], 
                   Row[{"n=", n, 
                     Style[Row[{"(", N[100 n/Nn, {2, 9}], "%)"}], Small, Italic, Gray]}]}],
               Epilog -> {{Red, PointSize[Medium], Point[points]}}],

           (* otherwise display a message *)
           Column[{Row[{"No points"}], Row[{"f=", N[pL, {2, 9}]}], Row[{"r=", N[r, {2, 9}]}]}]]]]]]
  ),
  SynchronousUpdating -> False,
  SynchronousInitialization -> False,
  Deployed -> True
]

The legend of the plot displays the current values of the controls (f and r) and the number of outliers (n) currently identified (both in absolute numbers and as a percentage of the total number of elements in data).

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.