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

Question

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.

Thanks in advance for your help! Tom

Answer 1

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

Answer 2

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 *)

Answer 3

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 *)

Answer 4

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 = Developer`PartitionMap},
    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

Answer 5

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.