Values of counting functions

Question

I have a table of integers, say A, and consider the counting function C, which, for a given x, gives the number of elements in A which are less or equal to x. Next, I want to plot C(x) in some range, say for x less then or equal to N. The question is: How this can be done quickly?

A very naive approach is illustrated by the following example.

max = 10^3;
n = 10^3;
A = RandomInteger[{1, max}, n];
c[k_]:=c[k]=Length[Intersection[A,Range[k]];

and c[k] gives the number of integers in A which are less then or equal to k. One can also use

cc[k_]:=cc[k]=Count[A, u_ /; u < k + 1];

However, if the set A has, say n = 10^5 elements, then to compute all values of c[k] is time consuming. I am pretty sure that there should be a quick and efficient method to compute all values of c[k] in a given range, but I failed to find one.

Answer 1

Here is another way to go if the list is only nonnegative integers.

lessEqualCount[list_, max_] := Module[
  {st, acc = 0, j = 1},
  st = SortBy[Tally[list], First];
  Table[If[j <= Length[st] && st[[j, 1]] == k, acc += st[[j, 2]]; j++];
   acc,
   {k, max}]
  ]

I'll show on a small example.

In[24]:= SeedRandom[1234];
ll = RandomInteger[{1, 10}, 20]

(* Out[25]= {1, 7, 10, 7, 1, 8, 1, 1, 9, 5, 5, 9, 6, 10, 8, 3, 9, 5, 6, 9} *)

There are four 1's, no 2's, 1 3, no fours, three 5's,.... So the result should begin {4,4,5,5,8,...}.

In[26]:= lessEqualCount[ll, 10]

(* Out[26]= {4, 4, 5, 5, 8, 10, 12, 14, 18, 20} *)

This gives an idea of speed.

max = 10^5;
n = 10^6;
alist = RandomInteger[{1, max}, n];
Timing[lessEqualCount[alist, max];]

(* Out[30]= {0.20433, Null} *)

Possibly this could be sped up by using Compile. Alternatively one might make a 0-order interpolation, using only the values in the tally and skipping the ones that don't appear in the list.

--- edit ---

This variant is a bit better in that it allows to use a minimum start not equal to 1 and allows to have values in the list less than that minimum.

lessEqualCount2[list_, min_, max_] := Module[
  {st, acc = 0, j = 1},
  st = SortBy[Tally[list], First];
  While[j <= Length[st] && st[[j, 1]] < min, acc += st[[j, 2]]; j++];
  Table[If[j <= Length[st] && st[[j, 1]] == k, acc += st[[j, 2]]; j++];
   acc,
   {k, min, max}]
  ]

--- end edit ---

Answer 2

You can use SparseArray and Accumulate to assemble a vector cvec so that cvec[[k]] == c[k]. When A has length 10^5, then creating cvec is already faster than a single call to Length[Intersection[A,Range[max]]].

min = -10^3;
max = 10^3;
n = 10^5;
A = RandomInteger[{min, max}, n];


With[{
   (*Remember the current SparseArrayOptions.*)
   spopt = SystemOptions["SparseArrayOptions"]
   }, 
   Internal`WithLocalSettings[
     (*Modify SparseArrayOptions so that SparseArray uses additive assembly.*)
     SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> Total}];,(*Do the actual assembly.*)
     v = SparseArray[Partition[A - min + 1, 1] -> 1, {max - min + 1},0];
     ,
     (*Set SparseArrayOptions to previous value.*)
     SetSystemOptions[spopt]];
   ];

(*Convert to a dense vector and accumulate.*)
cvec = Accumulate[Normal[v]]; // AbsoluteTiming

This does the same as the following code, only 30 times faster:

v = ConstantArray[0, max - min + 1];
Do[++v[[a-min + 1]], {a, A}];
cvec = Accumulate[v];

Caveat: This works only if A does not contain any integers $< 1$. You might have to add an offset to A otherwise.

Edit

Here is a compiled and parallelized version of the above Do-loop approach. Since SparseArray is not parallelized and does some extra work, this will work approximately $10 \times \mathrm{no of threads}$ faster.

cAccumulatedCounts = 
  Compile[{{A, _Integer, 1}, {min, _Integer}, {max, _Integer}, {threadcount, _Integer}, {thread, _Integer}},
   Block[{offset, counts, n, a, b, begin, end},
    n = Length[A];
    begin = 
     Quotient[n, threadcount] (thread - 1) + 
      Mod[n, threadcount] Quotient[thread - 1, threadcount] + 1;
    end = 
     Quotient[n, threadcount] thread + 
      Mod[n, threadcount] Quotient[thread, threadcount];
    
    offset = 1 - min;
    counts = Table[0, {max + offset}];
    
    Do[++counts[[Compile`GetElement[A, i] + offset]], {i, begin, end}];
    Accumulate[counts]
    ],
   CompilationTarget -> "C",
   RuntimeAttributes -> Listable,
   Parallelization -> True,
   RuntimeOptions -> "Speed"
   ];

AccumulatedCounts[A_?(VectorQ[#, IntegerQ] &), min_Integer, max_Integer] := Module[{threadcount},
   threadcount = "ParallelThreadNumber" /. ("ParallelOptions" /. SystemOptions["ParallelOptions"]);
   Total[
    cAccumulatedCounts[A, min, max, threadcount, Range[threadcount]]]
   ];

AccumulatedCounts[A_?(VectorQ[#, IntegerQ] &)] :=  AccumulatedCounts[A, Sequence @@ MinMax[A]];

If you happen to know the bounds of the array entries, then you can do

c = AccumulatedCounts[A, min, max];

Otherwise call the following; this will compute min and max first:

c = AccumulatedCounts[A];

Answer 3

Here is an alternative method which should work well even if negative numbers are present.

Let's generate some input data. Note this is about 8GB of data just for the starting list, so if you have less than probably 25GB of RAM available then you may need to lower n for this methodology.

n = 10^9;
alist = RandomInteger[{-1000,1000}, n];

The idea here is that the counting function of "how many elements are less than or equal to x in the list?" is equivalent to "what is the last position containing x in the sorted list?".

To find that efficiently we count each instance of x using Tally[alist], and then sort this tally (which is a list of {x, n} pairs indicating that we found n copies of x) by the x's.

SortBy[Tally[alist], First]

Using #[[1]]& to sort by only the x values. SortBy will arrange the list based on whatever this function returns in ascending order, and this function just returns the first value of each pair.

We can do a running total of the elements of this Tally using Accumulate, but Accumulate applies to all columns by default. Thus, we'll transpose the data twice and accumulate in the middle:

Transpose[{First[#],Last[#]}&[Transpose[
    SortBy[Tally[alist], First]
  ]]

This gets us a list of {x, c[x]} values. We can put a lower limit on this list by Joining a k such that c[k] is 0 to the front:

Join[{{Min[alist] - 1, 0}},
    Transpose[{First[#],Last[#]}&[Transpose[
        SortBy[Tally[alist], First]
      ]]]

And this we can put into an Interpolation of order 0 to get a functional form. A zeroth order interpolation simply fits everything to the immediately succeeding point in the fitting curve, creating a perfectly level step function which goes the wrong direction with the step to match c[x].

This can be dealt with for fairly minimal performance loss though, we'll just have to add 2 manual data points to the ends of the data set and reverse the x values twice:

c = Evaluate[
      Interpolation[
        Join[{{-1 - Max[alist], Length[alist]}}, 
             Transpose[{-First[#], Accumulate[Last[#]]} &[
               Transpose[SortBy[Tally[alist], First]]]],
             {{1 - Min[alist], 0}}],
        InterpolationOrder -> 0][-#]] &

Wrap the above in AbsoluteTiming, start from a fresh kernel for timing, and we get:

n = 10^9;
alist = RandomInteger[{-1000, 1000}, n];
AbsoluteTiming[c = Evaluate[
      Interpolation[
        Join[{{-1 - Max[alist], Length[alist]}}, 
             Transpose[{-First[#], Accumulate[Last[#]]} &[
               Transpose[SortBy[Tally[alist], First]]]],
             {{1 - Min[alist], 0}}],
        InterpolationOrder -> 0][-#]] &]

{ 2.98903, InterpolatingFunction[...] }

Seems decently efficient to me.

c[0] should be around 500,000,000 depending on the data set you specifically generated. If you put SeedRandom[17] in front, c[0] returns 500250556.

Answer 4

accumulatedBinCounts = Rest @ Accumulate @ BinCounts[#, 1] &;

Using Daniel's input examples:

SeedRandom[1234];
ll = RandomInteger[{1, 10}, 20];

accumulatedBinCounts @ ll

{4, 4, 5, 5, 8, 10, 12, 14, 18, 20}

max = 10^5;
n = 10^6;
SeedRandom[1];
alist = RandomInteger[{1, max}, n];


First @ Timing[res1 = lessEqualCount[alist, max];]

0.831814

First @ Timing[res2 = accumulatedBinCounts @ alist;]

0.042896

res1 == res2

True

Answer 5

You can use Sort with Split. This method works for both positive and negative integers.

With

max = 10^5;
n = 10^5;
SeedRandom[n];
a = RandomInteger[{1, max}, n];

Then pairs containing the integer and its count can be obtained by

res = Map[Through[{First, Length}@#] &, Split@Sort@a];
Short@res

{{2,2},{3,1},{6,1},<<63042>>,{99997,1},{99998,1},{100000,2}}

where the first item is the integer and the second its count.

There are a few ways to access the data.

Select can give you the pair if it exist and an empty list when it does not.

Select[First@# == 6 &]@res

{{6, 1}}

Select[First@# == 4 &]@res

{}

Converting res to an Association will give access by key (the integer) but return Missing for keys that do not exist.

ares = Association[Rule @@@ res];
ares[6]

1

ares[4]

Missing["KeyAbsent", 4]

However, you can merge ares with an Association of defaults to return a value on all integers.

allres = Merge[
   {
    AssociationThread[Range@n, 0]
    , Association[Rule @@@ res]
    }
   , Total
   ];

Now

allres[6]

1

allres[4]

0

Hope this helps.