Skip to main content

probability or statistics - How to remove outliers from data


I have a data with noise which some times includes significant outliers. The position of the outliers are random.


For example:


data1 = Table[PDF[NormalDistribution[3.5, .8], i], {i, -5, 15, .01}] +
RandomReal[{100, 500}];
noise = RandomReal /@ RandomReal[{-0.2, .2}, Length[data1]];
data2 = data1 + noise;
n = RandomInteger[{1, Length[data2]}, RandomInteger[{2, 10}]];
data2[[n]] = data2[[n]]*1.01;
ListPlot[{data2}, PlotRange -> All]


enter image description here


One solution is to use the average of the data but because the position of the outliers are random, the non-noisy data is hard to extract. The whole level of the data is random which means I can not use fixed reference to check and remove outliers.


Any idea how to remove these points using Mathematica?


Thanks.



Answer



I will give you two similar methods. But, I will rewrite one of the comments above just to make sure it is read.




You've been given some fine answers, but be absolutely sure that removing the outliers is Doing The Right Thing™. You might want to consider "robust" methods that can deal with the presence of outliers. – Guess who it is.





Simple Gaussian Threshold


The simplest way is to remove the moving mean of the data, then compute its standard deviation (σ), then pick a level at which you want to reject the data, say at 1%, so you can remove any points that vary more than 3×σ . If you know how the data is distributed about its mean values, then you can pick a different method. You can also remove the median since that would be less sensitive to the distribution.


SeedRandom[1245]; 
data1 =
Table[PDF[NormalDistribution[3.5, .8], i], {i, -5, 15, .01}] +
RandomReal[{100, 500}];
noise = RandomReal /@ RandomReal[{-0.2, .2}, Length[data1]];
data2 = data1 + noise;

n = RandomInteger[{1, Length[data2]}, RandomInteger[{2, 10}]];
data2[[n]] = data2[[n]]*1.01;
ListPlot[{data2}, PlotRange -> All]

enter image description here


We have about 8 outliers. We compute the moving average,


movingAvg=ArrayPad[MovingAverage[data2, 5],{5-1,0},"Fixed"]

Here we subtract the moving mean,


subtractedmean = (Subtract @@@Transpose[{data2, movingAvg}]);


enter image description here


Now find the locations of the outliers:


 outpos=Position[subtractedmean, x_ /; x>StandardDeviation[subtractedmean]*3];

Length[outpos]


8




looks like we got the right number of outliers. Removing them.


 newdata=Delete[data2,outpos]
ListPlot[newdata, PlotRange -> All]

enter image description here


To give you an idea of "Threshold" line in this case,


 dathreshold = 
ConstantArray[StandardDeviation[subtractedmean]*3, Length[data2]] + movingAvg;

Here is the "Threshold" line drawn along with the points removed,



 Show[ListPlot[data2, PlotRange -> All, AspectRatio -> 1], 
ListPlot[dathreshold, Joined -> True, PlotStyle -> {Thick, Purple}],
Graphics[{Red,Circle[#, {100, 0.5}] & /@
Thread[{First /@ outpos, data2[[First /@ outpos]]}]}]]

enter image description here


By derivatives


A second way to remove outliers, is by looking at the Derivatives, then threshold on them. Differences in the data are more likely to behave gaussian then the actual distributions.


 diff=Abs@Differences[data2,2];


ListPlot[diff, PlotRange -> All, Joined -> True]

enter image description here


Now you do the same threshold, (based on the standard deviation) on these peaks. Note that the outliers are now really well separated from the actual data. You can find the peak positions that are above the threshold you set, in our case we will keep using 3×σ.


You can probably use the peak finding function from V10 (not sure if there is a way to threshold the peaks), but since I stuck in V9 I do the poor's man way.


 newpos=Flatten[Position[Partition[diff, 3, 1], 
x_ /; ((x[[1]] < x[[2]] > x[[3]]) && (x[[2]] > 3*sddif)), {1}]] + 2

newpos===First /@ outpos



True



Comments

Popular posts from this blog

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...