Fixing the problem with HistogramList boundary values

Question

I am using HistogramList to bin the data and get frequency counts to generate probability distributions. I've noticed something surprising and terrible about Mathematica 10.4's binning method: it leaves out any points on the highest bin border.

TestData={1,1,1,1,2,2,2,3,3,3,3,4,4,4,5,5,5,6,6,7,7,8,8,9,9};
NumberOfBins=10;
TheMinMax=MinMax[TestData];
HistogramBinWidth=N@(TheMinMax[[2]]-TheMinMax[[1]])/NumberOfBins;
HistogramData=HistogramList[TestData,{TheMinMax[[1]],TheMinMax[[2]],HistogramBinWidth}]

The results of this is:

{{1.,1.8,2.6,3.4,4.2,5.,5.8,6.6,7.4,8.2,9.},{4,3,4,3,0,3,2,2,2,0}}

You see that the counts for the value 9 are missing from the final bin. If you check:

Length[TestData]
Total[HistogramData[[2]]]

We see that there are 25 items in the list, but only 23 counts in the HistogramList output. This is clearly a mistake.

The reason is clearly that Mathematica uses a use bin boundaries like this: [a,b),[b,c),...which works up until the last bin. That needs a special handler to include those points too, and Wolfram forgot about it.

So the question is, what is the best way to compensate for this error/bug?

I've padded the min and max by a tiny amount to get it to work, but that introduces error (it's a kludge rather than a fix). One option would be to test the list for points that match the upper bound and manually add them to the counts for that bin. Another would be to bypass the HistogramList function completely and do the correct binning manually. But which of these is the fastest? What's the best way to get accurate bin counts?

Answer 1

As I mentioned above, one obvious solution is to test the list for points that match the upper bound and manually add them to the counts for that bin. Here is some basic code that does that:

AddToFinalBin=If[TheMinMax[[2]]==HistogramData[[1]][[-1]],Count[TestData,x_/;x==HistogramData[[1]][[-1]]],0]
HistogramData[[2]][[-1]]=HistogramData[[2]][[-1]]+AddToFinalBin

Note that the conditional is there so that the same bins can be used for different lists with values within the same range, but not necessarily reaching the max (i.e., boundaries from a whole time series, but binned separately by time window).

Perhaps this code can be improved, and especially this is easy in one dimension, but seems like an inefficient approach for 2D data. Perhaps there is an entirely different and better way.

Answer 2

Here's a version of the other solution, doing the binning and counting from "scratch".

HistogramBinWidth=(TheMinMax[[2]]-TheMinMax[[1]])/NumberOfBins;
TransitionsPairs=Partition[Riffle[Drop[TheData,-1],Drop[TheData,1]],2];    
(* this is just an easy way to convert the 1D series data into 2D data *)

BinBoundaries=Table[TheMinMax[[1]]+(i*HistogramBinWidth),{i,0,NumberOfBins}]; 
(* these are the bins for one dimension, but they are the same in both dimensions *)

FrequencyData=ConstantArray[0,{NumberOfBins,NumberOfBins}];
Do[
    If[BinBoundaries[[i]]<=k[[1]]<BinBoundaries[[i+1]]&&BinBoundaries[[j]]<=k[[2]]<BinBoundaries[[j+1]],FrequencyData[[i]][[j]]++,##&[]];
    If[i==NumberOfBins&&k[[1]]==BinBoundaries[[i+1]]&&BinBoundaries[[j]]<=k[[2]]<BinBoundaries[[j+1]],FrequencyData[[i]][[j]]++,##&[]];
    If[j==NumberOfBins&&BinBoundaries[[i]]<=k[[1]]<BinBoundaries[[i+1]]&&k[[2]]==BinBoundaries[[j+1]],FrequencyData[[i]][[j]]++,##&[]];
    If[i==NumberOfBins&&j==NumberOfBins&&k[[1]]==BinBoundaries[[i+1]]&&k[[2]]==BinBoundaries[[j+1]],FrequencyData[[i]][[j]]++,##&[]];
,{i,NumberOfBins},{j,NumberOfBins},{k,TransitionsPairs}]; 

This basically just handles the cases logically separately, but done together in the same loop. I tried it with nested If statements and the Timing was the same. This is very much a brute-force approach, so although it works fine I expect there is a more elegant way to do the same thing.

Answer 3

Using the method of building from scratch totally works in terms of getting the correct frequency data, but it is MUCH slower than the built-in histogram function: Timing = 0.0312 vs Timing = 0..

Since the bug only excludes values exactly equal to the highest bin boundary, an intermediate solution is to specify the bins first, then add padding to the highest bin boundary, then use the built-in function.

BinBoundaries=Table[TheMinMax[[1]]+(i*HistogramBinWidth),{i,0,NumberOfBins}];
BinBoundaries[[-1]]=TheBins[[-1]]+.0000001;
HistogramData=HistogramList[TransitionsPairs,{BinBoundaries}]

I used a tiny amount of padding, but it doesn't matter for me because this only extends the highest bin, which (in my case) is set to the largest value in the data. So in my case it can't erroneously include points outside the desired range.

Of course, if in other applications the upper bound is not the max datapoint then this solution may not work properly, and Wolfram still needs to fix the bug.