Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question

5 Answers 5

up vote 8 down vote accepted

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

share|improve this answer
    
Your 2 samples give different answer… –  xzczd Sep 4 '13 at 2:35
2  
In fact you can also use the 3rd argument of Select to take the first matched sublist. –  xzczd Sep 4 '13 at 3:26
1  
@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? –  Mr.Wizard Sep 4 '13 at 15:22
1  
@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. –  Jonie Sep 5 '13 at 11:38
1  
@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. –  Mr.Wizard Sep 5 '13 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 *)
share|improve this answer

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 *)
share|improve this answer
    
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. –  Mr.Wizard Sep 4 '13 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 = 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
share|improve this answer
    
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). –  Michael E2 Sep 4 '13 at 11:43
    
+1. I always learn something new from your answers. –  RunnyKine Sep 5 '13 at 12:48
    
@RunnyKine Thanks; I try! (A bit too hard in this case perhaps; I think my verbosity put off Tom.) –  Mr.Wizard Sep 5 '13 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.

share|improve this answer
    
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. –  Mr.Wizard Sep 4 '13 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. –  jVincent Sep 4 '13 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. –  jVincent Sep 4 '13 at 8:36
    
Select has the short-circuit behavior. That's the nice thing about it: you can easily go either way. –  Mr.Wizard Sep 4 '13 at 8:40

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.