6
$\begingroup$

I'm just trying to read in some data before I do some calculations. The data set is over 200,000 values that look like this:

Blockquote

Here is the code I am using:

ClearAll["Global`*"];

Hrank = Flatten[Import["C:\\Projects\\Points_Analysis\\H_Rank.txt", "Table"]];
HighRank = Cases[Hrank, x_?NumericQ /; x <= Mean[Hrank] + 3*StandardDeviation[Hrank] && 
 x > Mean[Hrank] - 3*StandardDeviation[Hrank]];

It's just determining the values that lie within 3 sigma of the mean, but the calculation is taking hours (so far 6!) and doesn't seem to want to end. Could someone please tell me what is going wrong?

$\endgroup$
16
$\begingroup$

Your code is recalculating Mean[Hrank] and StandardDeviation[Hrank] with each comparison, making it extremely slow. Store them separately, and the calculation becomes much faster. I simulated an Hrank vector using:

Hrank = RandomInteger[{0, 1000}, 200000];

And then did:

mean = Mean[Hrank];
stddev = StandardDeviation[Hrank];
HighRank = 
   Cases[Hrank, 
    x_?NumericQ /; 
     x <= mean + 3*stddev && x > mean - 3*stddev]; // AbsoluteTiming

All as one cell, so that it would automatically update if Hrank changed. This gave a timing (via AbsoluteTiming) of 3.421 seconds (on my system) for 200,000 elements. Using a smaller sample, it appears that this should give exactly the same results as the original code, as at no point would Hrank be altered during this.

It's possible to go faster than this using Compile and some cleverness. Since it's actually more elegant than using Cases, it's presented below:

Define a family of test functions (testf[m,s]) below:

testf[m_, s_] := 
  Compile[{{x, _Real}}, x <= m + 3 s && x > m - 3 s, 
   RuntimeAttributes -> {Listable}, Parallelization -> True];

Note that testf[m,s] returns a compiled function object which tests a real number as to whether it is within 3 standard deviations (True) or not (False). This is then used with Pick to pick out the final set:

HighRank2 = Pick[Hrank, testf[mean, stddev][Hrank]]; // AbsoluteTiming

It is roughly 200 times faster (on my system, probably faster still if you have more than 8 threads available) than the above version using Cases.

@Shadowray provides an even faster version utilizing UnitStep, RealAbs, and numericizing the input list to machine numbers before hand. Some cursory examination suggests the results match identically in at least most cases, so I'll preserve it here below:

threeSigmaFilter[list_] := 
 With[{nlist = N[list]}, 
  Pick[list, 
   UnitStep[
    RealAbs[nlist - Mean[nlist]] - 3 StandardDeviation[nlist]], 0]]

This centers the whole data list on the mean, takes the absolute value of that, and subtracts 3 standard deviations. If this result is less than 0, then UnitStep returns 0, and Pick's 3rd argument tells Pick to keep it from the original list. It doesn't seem like Compile can be used to further speed this up. This isn't very surprising, UnitStep and RealAbs (or Abs, not much performance difference there it seems) are already Listable and quite fast as built-ins, and there is a time overhead to compiling a function.

$\endgroup$
  • 2
    $\begingroup$ Can be further optimized to something like: threeSigmaFilter[list_] := With[{nlist=N[list]}, Pick[list, UnitStep[RealAbs[nlist-Mean[nlist]] - 3 StandardDeviation[nlist]], 0] ] $\endgroup$ – Shadowray Jul 26 at 18:33
  • $\begingroup$ @Shadowray Thanks for pointing that out. I was hesitant to numericize the list because it introduces a little bit of inaccuracy, but it is fair to say that it's quite a bit faster. I'm always surprised at how fast UnitStep is too. $\endgroup$ – eyorble Jul 26 at 20:26
  • $\begingroup$ I know we're not supposed to express thanks, but I just wanted to let both of you know how much I appreciated your solutions! $\endgroup$ – Angus Jul 27 at 12:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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