# How to find first list element that differs from average of N previous elements by more than a given amount?

Given a list of real numbers, I need to write a function that returns the value of the first element in the list that differs from the mean of the previous $N$ elements by greater than $P$ standard deviations of those same elements. (I use mean and standard deviation as sample functions, I may decide to use something else.)

Clearly, I can write a loop to do this, but I'd like to use functional programming techniques to do this, and continue to learn better Mathematica programming style. Select[ ] looked promising, but I see no way to involve the previous list elements in the selection criterion.

Here's an attempt (without select)

lst = {1, 2, 2, 1, 2, 5, 2, 4};
n = 4;
p = 1;
Flatten[If[Abs[Take[#, -1] - Mean[Drop[#, -1]]][[1]] > p *StandardDeviation[Drop[#, -1]],
Take[#, -1], {}] & /@ Partition[lst, n + 1, 1, 1]] // First


Alternatively with Select.

Select[Partition[lst, n + 1, 1, 1],
Abs[Take[#, -1] - Mean[Drop[#, -1]]][[1]] > p*StandardDeviation[Drop[#, -1]] &]
// First // Last


For a more efficient implementation - see Mr Wizard's answer

• Your 2 samples give different answer… Sep 4, 2013 at 2:35
• In fact you can also use the 3rd argument of Select to take the first matched sublist. Sep 4, 2013 at 3:26
• @Tom with my sincere apologies to Jonie, this is really not the best method to use. Far better would at least be PartitionMap so as not to waste massive amounts of memory. Have you seen the performance metrics in my answer? Have you taken the time to run your own? Sep 4, 2013 at 15:22
• @TomDickens, if you could kindly change the acceptance to Mr Wizard's answer for future visitors, that'd be great. Mr Wizard - no worries. I did learn that partition construct off another answer you made afterall. Sep 5, 2013 at 11:38
• @Jonie Thanks, I appreciate that. Partition certainly has many good uses. I'll add a new section to my (already too long, apparently) answer showing how you can use PartitionMap efficiently. Sep 5, 2013 at 11:48

Here is a functional programming style implementation based on Fold, using Throw and Catch to mimic the non-existent FoldWhile construct. You might be interested in some earlier attempts to build a FoldWhile function.

firstOutlier[data : {__?NumericQ}, p_?Positive, k_Integer?Positive] /;
k < Length[data] :=
Catch[Fold[
If[Abs[#2 - Mean[#1]] > p*StandardDeviation[#1], Throw[#2],
Join[Rest[#1], {#2}]] &, Take[data, k], Drop[data, k]] ]


Using Jonie's test data:

 firstOutlier[{1, 2, 2, 1, 2, 5, 2, 4}, 1, 4]

(* 5 *)


Here's a way using NestWhile:

f[lst_?(VectorQ[#, NumericQ] &), n_Integer?Positive, p_?Positive] :=
If[Length @ # > n, #[[n + 1]], None] &@
NestWhile[Rest, lst,
Length @ # > n &&
Abs[Mean @ #[[;; n]] - #[[n + 1]]] <= p StandardDeviation @ #[[;; n]] &]

f[{1, 3, 2, 0, 2, 2, 1, 2}, 4, 1]
(* None *)

f[{1, 3, 2, 0, 2, 5, 1, 2}, 4, 1]
(* 5 *)

• I wrote my answer before reading any of the others but the body of this certainly looks similar to my test. I guess that's the natural way to approach this however. Now for timings. Sep 4, 2013 at 7:20

### Considerations

This question reminds me of Iterate until condition is met in a way. You would probably find my answer there applicable as you will likely face similar considerations, e.g. block-based processing comes with a certain overhead but allows for vector optimizations in many cases.

Things like Mean will be more efficiently calculated in an incremental fashion using e.g. MovingAverage or ListCorreclate. Therefore there is a trade-off between pre-processing the data in a faster way and processing only as much as is needed until the first result is found, with the average length of your data and the position of the first match determining which is best. I don't know how to (efficiently at least) calculate an incremental standard deviation but that is likely my own ignorance.

## Basic proposal

For generic functions without the incremental approach I propose using Select on list of indexes with a test function that extracts elements from the input list and performs the requested comparison. This has considerable memory advantage over using Partition and will also be faster when the match occurs early in the input list.

(One could use an incrementing position variable rather than a list of indexes, which would conserve memory further. I chose to use the indexes because of the flexibility it provides such as easily returning the first m matches.)

Here is a partially generalized function that implements this proposal:

find[list_, n_, p_, m_: 1, f1_: Mean, f2_: StandardDeviation] :=
Module[{test},
test[x_] := Abs[list[[x]] - f1@#] > p * f2@# & @ list[[x - n ;; x - 1]];
list[[ Select[Range[n + 1, Length@list], test, m] ]]
]

• The three required parameters are list, n, and p.
• m is the maximum number of matches to return, defaulting to one.
• f1 and f2 default to Mean and StandardDeviation but could be replaced with any other vector functions.

To make find more general it would probably be best to make p part of the f2 function but I wanted to respect your defined parameters.

Example:

SeedRandom[1]
list = RandomReal[{-9, 9}, 50];

find[list, 5, 1.5]

{3.60853}

find[list, 5, 1.5, 3]

{3.60853, 4.47582, 8.58909}


### Tuning

For maximum speed but less flexibility one can Compile a simple Do loop. It is several times faster and has the same computational complexity as the method above.

Limitations:

• You lose the on-the-fly ability to change functions (f1, f2), at least without meta-programming to generate a new compiled function.

• It returns only the first match rather than a specified number (m), though that could be added with more code and e.g. Sow and Reap.

• It will only work on lists of machine-size numbers.

Code:

findC =
With[{f1 = Mean, f2 = StandardDeviation},
Compile[{{list, _Real, 1}, {n, _Integer}, {p, _Real}},
Do[
If[#, Return @ list[[x]]] &[
Abs[list[[x]] - f1@#] > p*f2@# &@list[[x - n ;; x - 1]]
],
{x, n + 1, Length@list} ]]];


### Alternative using PartitionMap

I favor my find because of the ease of specifying the number of matches to return while still preserving the short-circuit behavior, however it may be easier to think of this problem in terms of Partition.

But there are two major drawbacks to using Partition:

• you will expand the array to $N$ times its original size
• the entire vector is partitioned in advance of testing, causing a fixed overhead even on very early matches

Both these problems can be solved by using PartitionMap and Return, as I shall illustrate. Note that I will be using a second argument in Return.

With the test function from Jonie's Select method, refined, and some data:

test = Abs[Last[#] - Mean[Most @ #]][[1]] > p*StandardDeviation[Most @ #] &;

SeedRandom[1]
lst = RandomReal[{-9, 9}, 150000];
{n, p} = {20, 3.0};


This is equivalent to Jonie's second method (with the refined test):

Select[Partition[lst, n + 1, 1, 1], test] // First // Last // Timing

{0.827, -8.73536}


This match occurs only 864 places into the list; the slow speed is due partly to the Partition overhead but primarily to the fact that all places were tested, because the third parameter of Select was not used. (The memory consumption would remain high either way.)

Now here is the improved method:

Needs["Developer"]

PartitionMap[If[test@#, Return[#, PartitionMap]] &, lst, n + 1, 1, 1] // Last // Timing

{0.005488, -8.73536}


This is faster than the flexible find but not as fast as findC (timeAvg code below):

find[lst, n, p]  // timeAvg
findC[lst, n, p] // timeAvg

0.01436

0.001872


Finally, as a function:

findP[list_, n_, p_, f1_: Mean, f2_: StandardDeviation] :=
Module[{toss, PMap = DeveloperPartitionMap},
toss = If[Abs[Last[#] - f1[Most@#]] > p*f2[Most@#], Return[#, PMap]] &;
PMap[toss, list, n + 1, 1, 1] // Last
]


## Timings

I wrote this answer before looking at others so as not to be influenced. I see now that Michael used code very similar to my test function but he did not an indexed-based Select, choosing instead repeated shortening of the list with Rest. The Rest approach will slow down proportionally to the length of the input list because of the reallocation that occurs with Mathematica lists while the index approach will not. Jonie's code pre-partitions the entire list an will use n times the memory of the input list so I did not attempt to use it on the long lists below.

Correction: Contrary to what I stated before, Verbeia's code does not suffer from the reallocation problem.

On short lists the functions are in the same magnitude, but on long lists:

SeedRandom[1]
list = RandomReal[{-9, 9}, 500000];

f[list, 5, 9]    // Timing  (* Michael's function *)
find[list, 5, 9] // Timing

{2.668, -4.43862}

{0.047, {-4.43862}}


The performance gap increases as the list gets longer:

SeedRandom[1]
list = RandomReal[{-9, 9}, 1500000];

f[list, 5, 20]    // Timing
find[list, 5, 20] // Timing

{46.379, -4.75616}

{0.203, {-4.75616}}


Now fincC and findP on the same data:

findC[list, 5, 20] // Timing

findP[list, 5, 20] // Timing

{0.058, -4.75616}

{0.156, -4.75616}


### Timings, part 2

Correcting myself, Verbeia's code does not suffer from the reallocation problem described above. Here are some comparative timings against my code (timeAvg posted many times before).

On a test that results in early termination find is and order of magnitude faster than firstOutlier, and findC is another three orders faster than that:

SetAttributes[timeAvg, HoldFirst]

timeAvg[func_] := Do[If[# > 0.3, Return[#/5^i]] & @@ Timing@Do[func, {5^i}], {i, 0, 15}]

{size, n, p} = {1000000, 100, 1.7};

SeedRandom[1];
list = RandomReal[{-9, 9}, size];

firstOutlier[list, p, n] // timeAvg
find[list, n, p]         // timeAvg
findP[list, n, p]        // timeAvg
findC[list, n, p]        // timeAvg

0.327

0.02496

0.003864

0.000023936


However, with late termination find falls behind firstOutlier.

{size, n, p} = {1000000, 100, 2.4};

SeedRandom[1];
list = RandomReal[{-9, 9}, size];

firstOutlier[list, p, n] // Timing
find[list, n, p]         // Timing
findP[list, n, p]        // Timing
findC[list, n, p]        // Timing

{5.055, 8.71491}

{5.569, {8.71491}}

{4.929, 8.71491}

{2.09, 8.71491}


### Memory

Here is an illustration of the problem with Partition that I mentioned. Jonie's second method tested a bit better than the first so I will use that:

jonie2[lst_, n_, p_] :=
Select[Partition[lst, n + 1, 1, 1],
Abs[Take[#, -1] - Mean[Drop[#, -1]]][[1]] > p*StandardDeviation[Drop[#, -1]] &] //
First // Last


Now, each time in in a fresh kernel and starting with:

{size, n, p} = {500000, 100, 2.2};
SeedRandom[1];
list = RandomReal[{-9, 9}, size];


Mine:

find[list, n, p] // Timing
MaxMemoryUsed[]

{0.281, {-8.59958}}

37412144


Jonie's:

jonie2[list, n, p] // Timing
MaxMemoryUsed[]

{4.212, -8.59958}

923440192


findP corrects this, taking even less memory than find:

findP[list, n, p] // Timing
MaxMemoryUsed[]

{0.25, -8.59958}

27769328

• Thanks for the remark on Rest - I hadn't realized the memory issue. Your find beats my f hands down, but the timing complexity depends on the position of the first outlier. With a relatively small list of 150000, find was between 8.4 and 1.7 times faster (used SeedRandom[i] for i from 1 to 100). Sep 4, 2013 at 11:43
• +1. I always learn something new from your answers. Sep 5, 2013 at 12:48
• @RunnyKine Thanks; I try! (A bit too hard in this case perhaps; I think my verbosity put off Tom.) Sep 5, 2013 at 12:50

Fun little question. It seems to be getting a lot of good answers, but in any case I'll throw in my 2 cents. My immediate though is that this sort of fits the loop construct LengthWhile, since we really just want to keep on counting until we get a hit, but it doesn't naturally work on trailing lists, so we'll just add this to it. In that regard here's my suggesting:

countWhileList[list_, n_, criteria_] :=
LengthWhile[
Range[n + 1, Length@list],
criteria[list[[# - n ;; #]]] &
] + 1 + n


Which can then be used as:

normalDeviation[list_] := Abs[(Mean[list] - list[[-1]])/StandardDeviation[list]]
countWhileList[testData, 40, normalDeviation[#] < 1 &]


Which will count trailing lists of 40 elements and return the number of the first hit that fails the criteria. Depending on preference, it could instead be made to count the number of successes by removing the +1.

Naturally if you then want the value, you can use Part as in

Part[testData,countWhileList[testData, 40, normalDeviation[#] < 1 &]]


Keeping in this style but allowing for multiple returns you could use instead

selectWhereList[list_, n_, criteria_, m_] :=
Select[Range[n + 1, Length@list], criteria[list[[# - n ;; #]]] &, m]

Part[test, First@selectWhereList[test, 40, normalDeviation[#] > 3 &, 1]]


Which is essentially just the inner part of Mr. Wizards solution. Personally I prefer my first, if only because I so rarely get to use LengthWhile for anything constructive.

• This appears to be the same approach I took but rephrased to use LengthWhile instead of Select. I prefer Select because of the ability to return multiple matches. Sep 4, 2013 at 8:06
• @Mr.Wizard Thanks for the correction. My notion was that getting the value of an element in a list when you have the index is trivial, but it doesn't hurt to add. Sep 4, 2013 at 8:21
• @Mr.Wizard Also, while I understand a preference for Select in cases where you want all matches, it seems detrimental to this case which is specifically targeted at a first find, and therefore benefits from a short circuit and therefore a loop construct. Sep 4, 2013 at 8:36
• Select has the short-circuit behavior. That's the nice thing about it: you can easily go either way. Sep 4, 2013 at 8:40