Replacing a for-loop with something functional

I'm a physics grad and have been using Mathematica mostly for data visualization. The evaluation part I normally do with C/Python, but I've been wanting to try out Mathematica for these tasks, too.

Case in point, I made two notebooks for trimming and evaluating data for a recent experiment.

Code example:

new = data[[All, 1]]
For[i = 1, i < Length[data[[All, 1]]] - 1, i += 1,
If[
Abs[data[[i, 1]] - data[[i + 1, 1]]] > 2000,
new[[i + 1]] = Mean[{data[[i, 1]], data[[i + 2, 1]]}];
i += 1;
]];


data[] is a list with ~30000 datapoints (x, y, value). Here, x from xmin to xmax has to be linear, but there are some (very few) points that were measured incorrectly and thus are way above/below the other points.

The loop above iterates over this list and checks the difference of list element [i] and [i+1]. If this difference is larger than a threshold, element [i+1] is replaced by the mean value of element [i] and [i+2].

I would love to know if there was a more elegant way of doing this? Ideally faster, too, since I intend to use this on much larger lists. The above is one loop of many, so an explanation of the thought process that went into the solution would help a great deal :)

EDIT: Maybe I haven't stated my problem clearly enough, so I'll try again. I really don't need help regarding what the above code does - it does exactly what I want and has been tested with many datasets. I just want to know if there is a more "Mathematica-ish" way of writing it, since I've been told to avoid for-loops as Mathematica promotes a more functional style of programming.

ONE MORE EDIT: I'd like to thank everyone who devoted some time to this. I see that I should've included sample data when stating the problem, I'll try to remember next time ;)

• Why the -1 in Length[data[[All, 1]] - 1]? – enzotib Apr 1 '15 at 17:39
• Because of Mean[{x[i], x[i+2]}] - I'm cheating a bit because I'm not actually checking the last two elements of the list. Which is okay in this case, though. – mrbaozi Apr 1 '15 at 17:48
• You need Length[data[[All, 1]]] - 1, i.e. the -1 out of the Length call. – enzotib Apr 1 '15 at 17:59
• Woooooops yes you're right, damn brackets. Do need the -1 though, otherwise I'd have to do i <= length -2... – mrbaozi Apr 1 '15 at 18:05

Best to try and work with the entire lists rather than loop through it with For or Table.

You obviously have your reasons for your algorithm but to me it is strange because if your test returns True then you effectively increment twice. Which you can see with this version of your code:

list = {};
new = data[[All, 1]]
For[i = 1, i < Length[data[[All, 1]]] - 1, ++i,
test = Abs[data[[i, 1]] - data[[i + 1, 1]]] > 2000;
If[test,
list = Join[list, {{i, test}}];
new[[i + 1]] = Mean[{data[[i, 1]], data[[i + 2, 1]]}];
++i,
list = Join[list, {{i, test}}]]
];
list


So you are not actually looping through every element in your list.

Going through what you are doing

Abs[data[[i, 1]] - data[[i + 1, 1]]]


can be rewritten for the list as

Abs[-Differences[data[[All, 1]]]]


From here, If I wanted to "loop" through every element, rather than skip some, I would do this:

positions =
Pick[Range[2, Length[data[[All, 1]]] - 1],
Thread[Abs[-Differences[data[[;; -2, 1]]]] > 2000], True];

new[[positions]] =
Mean /@ Transpose[{data[[positions - 1, 1]],
data[[positions + 1, 1]]}];


So what happens in this code is that you test which elements need to be replaced in one go rather than via a loop. You then apply that test to find the positions in your list that get replaced. You could also use Position and I would recommend trying that as well to see which is the more efficient ...Pick is likely to be faster for longer lists. ListCorrelate is another function that could be used to develop a solution.

new1 = data[[All, 1]];
tmp = # > 2000 & /@ Abs[-Differences[data[[;; -2, 1]]]];
true = Rest@FoldList[If[#1 == #2, False, #2] &, False, tmp];
positions = Pick[Range[2, Length[data[[All, 1]]] - 1], true, True];
new1[[positions]] =
Mean /@ Transpose[{data[[positions - 1, 1]],
data[[positions + 1, 1]]}];
new1


Test

SeedRandom[123];
data = RandomReal[{1, 10000}, {30000, 2}];


On my computer -- V9.0.1 OS X -- your method (returning new) took 27 seconds mine (returning new1) took 0.054 seconds. A not too shabby 500 times improvement.

new==new1

(* True *)

• Thank you! The second method works. Would've loved to test the first one, but I get an error "symbol new does not have an immediate value". Not quite as short of a solution as I had expected, but it's definitely faster. I thought Mathematica had one-liners for everything ;). My code is a bit of a "hack", since ideally you'd want to fit the data and then compare each data point to that fit instead of a 2-point mean, but that isn't neccessary in this case. Also, I don't have a fit function for everything (data is only linear in this special case, the whole notebook has some 300+ lines of code). – mrbaozi Apr 2 '15 at 6:06
• @mrbaozi if you have 10.1 the SequencePosition approachlooks to be better ...but I have not tested it. – Mike Honeychurch Apr 2 '15 at 23:12

There are some problems with your example code, you are incrementing i twice, it looks like you might subscript past the end of data, etc.

Consider this idea. It sounds like you are looking at triples in your data. So start with

Partition[data, 3, 1]


That is going to take overlapping triples of your data and you process each triple separately.

Now write a function that takes three consecutive items and returns a (potentially) cleaned first item.

f[{{x1_,__},{x2_,__},{x3_,__}}]:=If[x1-x2>2000, (x1+x3)/2, x1];


That takes a list of three rows from your data and if the difference between the first item in the first two rows is too great then averages the first and third, otherwise it gives you the unchanged first item. (This doesn't exactly match your verbal description, but your verbal description doesn't match your example code, Hopefully you can see how you want to change this to get what you really want).

Now you would like to use that function f on every triple, so

Map[f, Partition[data, 3, 1]]


Given data = {{1875.81, 1, 2}, {-1613.51, 2, 3}, {533.669, 3, 4}, {-1758.47, 4, 5}, {1613.63, 5, 6}, {100.325, 6, 7}}

it give you

{1204.74, -1613.51, 1073.65, -1758.47}


and you can try to check if that does exactly what you want

• Thanks for this, looking into it! I can't guarantee Length[data] % 3 == 0, would that break the partitioning? Also, I need the resulting list to be of the same dimension as the original list, that's why I assign new = data[[All, 1]]. – mrbaozi Apr 1 '15 at 18:02
• @mrbaozi Look at the syntax of Partition carefully and check what the third argument does. – Szabolcs Apr 1 '15 at 18:31
• Okay, this is really fast, just tested it against my loop, which is cool. But it doesn't work. When the error is at index 2 in the partitioned list, it gets eliminated - intended behavior. In the next sublist, the error will be at index 1 though, and then the mean of the error and index 3 will be returned. This is why, in my loop, I increment the counter by 1 if an offset has been found. Any ideas? – mrbaozi Apr 1 '15 at 18:59
• Bill, your proposal seems to produce a wrong result. Example with data {{455, 1, 1}, {2453, 2, 2}, {2603, 3, 3}, {3102, 4, 4}, {2968, 5, 5}, {3373, 6, 6}, {139, 7, 7}, {2570, 8, 8}, {3432, 9, 9}, {801, 10, 10}} produces for @mrbaozi: {455, 2453, 2603, 3102, 2968, 3373, 5943/2, 2570, 3432, 801} and for Bill {455, 2453, 2603, 3102, 2968, 5943/2, 139, 2570}. You may crosscheck on this... – penguin77 Apr 2 '15 at 2:01
• @Szabolcs Partition creates a problem in the proposal of Bill, because the original code increments twice if the a value is replaced (once in the For-loop and once in the IF-statement. It means that the replaced value needs to be skipped.... – penguin77 Apr 2 '15 at 2:15

For demonstration purpose I'll use some simple test data

test = Range[20];
test[[5]] = 200;
test[[7]] = -200;

test


{1, 2, 3, 4, 200, 6, -200, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20}

Here is a way using one of the Sequence* functions introduced in version 10.1 of Mathematica:

(test[[# + 1]] = Mean@test[[{#, # + 2}]]) & /@
SequencePosition[test, {a_, b_} /; Abs[a - b] > 2000, Overlaps -> False][[All, 1]];

test


{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20}

new = data[[All, 1]];
(new[[# + 1]] = Mean@new[[{#, # + 2}]]) & /@
SequencePosition[new, {a_, b_} /; Abs[a - b] > 2000, Overlaps -> False][[All, 1]];


new = ReplacePart[new,
Thread[# + 1 -> Mean@new[[{#, # + 2}]] & /@ #] &
[SequencePosition[new, {a_, b_} /; Abs[a - b] > 2000, Overlaps -> False][[All, 1]]]]

• +1 for SequencePosition. I liked these new functions. – Murta Apr 1 '15 at 21:00
• The final method doesn't return the same output as his current procedural method, presumably because you are not capturing the skipped increments – Mike Honeychurch Apr 1 '15 at 23:49
• @ Karsten 7, I agree with Honeychurch and some problem as with Bill's problem – penguin77 Apr 2 '15 at 2:30
• @MikeHoneychurch I removed ReplaceRepeated approach. – Karsten 7. Apr 2 '15 at 9:41
• @penguin77 I removed ReplaceRepeated approach. – Karsten 7. Apr 2 '15 at 9:41
Table[
If[Abs[data[[i + 1]] - data[[i]]] > 2000,
Mean[{data[[i]], data[[i + 2]]}], data[[i]]],
{i, Length[data] - 2}]


It is unclear what you want to do with the last two data points, so the above code just eliminates them.

• Doesn't this include outliers twice? – mrbaozi Apr 1 '15 at 18:11