0
$\begingroup$

Here it is given a possible method in order to count data with weights. That solution works. However, it is ~140 times slower than the Bincounts function. See, e.g. the following code:

Ndata = 10^6;
SeedRandom[321];
data = RandomReal[{0, 10}, Ndata];
weights = RandomReal[{.9, 1.1}, Ndata];
bins = Table[i, {i, 0, 10}];

AbsoluteTiming[Last[HistogramList[WeightedData[data, weights], {bins}]]]
AbsoluteTiming[BinCounts[data, {bins}]]

I would like to know if it is possible to make a faster function.

A faster method could be to apply the function proposed here to the output of BinLists:

myFn = Merge[KeyIntersection[PositionIndex /@ {##}], Identity] &;
blist = BinLists[data, {bins}]
myFn[blist, Partition[data, 1]]
(*<|{0.986147} -> {{1}, {3}}, {1.49106} -> {{2}, {6}}, {3.23491} -> {{4}, {7}}, {7.15785} -> {{8}, {5}}, {8.9058} -> {{9}, {9}}|>*)

Then I could use the latter output (the positions) to sum the weights. However, myFn does not identify the positions if there are more elements inside a bin.

Does someone know how to improve upon this attempt? Other solutions are most welcome. Perhaps one could see how does the Python function numpy.bincount work?

$\endgroup$

1 Answer 1

3
$\begingroup$

I would use a solution based on Pick:

wBinCounts[data_, weights_, bins_] := Module[{min, max, b, d},
    {min, max} = MinMax[data];
    If[min<0 && max>0,
        d = data - min + 1;
        b = bins - min + 1,
        d = data;
        b = bins
    ];
    With[{ranges = Partition[b, 2, 1]},
        Total @ Pick[weights, Unitize@Clip[d, {##}, {0, 0}], 1]& @@@ ranges
    ]
]

For your data:

Ndata=10^6;
SeedRandom[321];
data=RandomReal[{0,10},Ndata];
weights=RandomReal[{.9,1.1},Ndata];
bins=Table[i,{i,0,10}];

AbsoluteTiming[Last[HistogramList[WeightedData[data,weights],{bins}]]]
AbsoluteTiming[BinCounts[data, {bins}]]
AbsoluteTiming[wBinCounts[data,weights,bins]]

{19.7562, {99667.3, 100028., 99756.9, 100022., 100024., 100361., 100513., 99787.1, 99992.4, 99881.1}}

{0.132353, {99688, 100024, 99735, 100031, 99998, 100371, 100488, 99778, 100001, 99886}}

{0.096202, {99667.3, 100028., 99756.9, 100022., 100024., 100361., 100513., 99787.1, 99992.4, 99881.1}}

Update

Clip[x, {min, max}] uses $min \leq x \leq max$, and so values of $x$ at the end points will be included in both intervals. A slightly slower version of wBinCounts that avoids this issue (using $min \leq x < max$):

wBinCountsLeft[data_, weights_, bins_] := With[{ranges=Partition[bins,2,1]},
    Total @ Pick[
        weights,
        BitXor[UnitStep[Subtract[data, #1]], UnitStep[Subtract[data, #2]]],
        1
    ]& @@@ ranges
]

And a version using $min <x\leq max$:

wBinCountsRight[data_, weights_, bins_] := With[{ranges=Partition[bins,2,1]},
    Total @ Pick[
        weights,
        BitXor[UnitStep[Subtract[#1,data]], UnitStep[Subtract[#2,data]]],
        1
    ]& @@@ ranges
]

Example:

wBinCountsLeft[{1,1.5,2}, {1,1,1}, {0,1,2,3}]
wBinCountsRight[{1,1.5,2}, {1,1,1}, {0,1,2,3}]

{0, 2, 1}

{1, 2, 0}

$\endgroup$
4
  • $\begingroup$ Thanks! Could you please generalize it to the case of bins that are not unit intervals? $\endgroup$
    – Valerio
    Jun 7, 2017 at 16:25
  • $\begingroup$ wow! it is even faster than BinCounts! $\endgroup$
    – Valerio
    Jun 7, 2017 at 22:38
  • $\begingroup$ It seems as if a data point is counted twice if it lies at the bin boundary. For example, wBinCounts[{1, 1.5}, {1, 1}, {0, 1, 2}] returns {1, 2} but there are only 2 points in total. Could you please fix this? $\endgroup$
    – Valerio
    Mar 16, 2018 at 13:20
  • $\begingroup$ thanks for the update! Could you please run the AbsoluteTiming of these new functions together with the old ones? I ran them on my Mac with Mathematica 11.3.0.0 and BinCounts got much faster. I believe I was using Mathematica 11.2.0.0 when I opened the question and I am not aware of any performance improvements regarding BinCounts. $\endgroup$
    – Valerio
    Mar 21, 2018 at 17:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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