np.histogram() inconsistent results

I think that #6100 may have introduced a bug in some edge cases. Here is a case found through fuzz testing:

|35> import numpy as np
|36> np.__version__
'1.10.4'
|37> arr = np.array([337, 404, 739, 806, 1007, 1811, 2012])
|39> hist, edges = np.histogram(arr, bins=8296, range=(2, 2280))
|40> mask = hist > 0
|41> left_edges = edges[:-1][mask]
|42> right_edges = edges[1:][mask]
|43> zip(arr, left_edges, right_edges)
[(337, 337.00000000000006, 337.27459016393448),
 (404, 404.00000000000006, 404.27459016393448),
 (739, 739.00000000000011, 739.27459016393448),
 (806, 806.00000000000011, 806.27459016393448),
 (1007, 1007.0000000000001, 1007.2745901639345),
 (1811, 1811.0000000000002, 1811.2745901639346),
 (2012, 2012.0000000000002, 2012.2745901639346)]

This was found through fuzz testing an internal accelerated histogram routine (which takes roughly the same strategy as #6100 but implemented in Cython) and comparing the results to np.histogram(). Comparing with numpy 1.9.2, everything worked identically. At least as of numpy 1.10.4, some of these edge cases popped up.

I believe the main difference between #6100 and our routine is that we do not precompute the histogram scaling factor bins / (mx - mn): https://github.com/numpy/numpy/blob/master/numpy/lib/function_base.py#L638
We compute ((tmp_a - mn) * bins) / (mx - mn) . The comments of our routine mention floating point problems as the reason for deliberately avoiding precomputation.

@astrofrog

Pinging @astrofrog

@rkern

@rkern Want to make a PR?

Hmm, looks like it might be linspace() to blame, actually.

Well, no, we use linspace() too. Hmm.

Ah, I see. We are doing some extra work to correct the location in case it happens to be off by one in these edge cases. Not sure if there is a fast way to do this vectorized in Python.

Sounds like I don't need to worry about this for 1.11.1 unless there is some obvious error involved.

The error is obvious: the general approach taken in #6100 has inherent inaccuracies. One could revert that PR, and the error will be fixed. The other fix (that mostly retains the performance gains of #6100) is obvious too: use the current method as just the first guess and adjust it up or down by 1 when needed. The only thing that isn't obvious is the precise implementation of the fix.

@rkern

@rkern - thanks for identifying this issue!

So just to be sure I understand, are you suggesting correcting the edges or the histogram values? I'm pretty sure it should be the former (the edges), but just want to make sure we are on the same page.

Conceptually, if the user requests:

np.histogram(arr, bins=8296, range=(2, 2280))

then what is being returned incorrectly is not the histogram but the edges, since in the first case, the edge position should be exactly 337 not 337.00000000000006, and if the user derives the exact solution, I think the histogram itself is correct currently. In other words, if we didn't return the edges, then the current histogram would be correct (and it's just that currently there is an inconsistency between the two).

It seems like linspace does have some floating point issues, since the following is incorrect:

In [3]: np.linspace(2, 2280, 8297)[1220]
Out[3]: 337.00000000000006

If we revert #6100, we will go back to a state where the histogram and edges are consistent but at the same time may be unintuitive due to the above issue (the bin edges won't be where the user expected them). If instead we fixed linspace, then I think histogram would then give the correct answers.

@mdickinson

The "histogram" is the combination of the bin counts and the bin edges. There is no histogram without the edges. Changing either such that they are consistent with each other fixes the problem; I don't think either approach is preferable to the other in terms of their result. If the user wanted fine control over the edges, they would have passed them in.

My suspicion is that you will not be able to change the bin edge computation in such a way that it is provably consistent with the indexing approach. I think we will have to still check for off-by-one errors and adjust the initial index if it points to a bin where the data point is not in its bounds. I'm happy to be proven wrong, though. @mdickinson, need an excuse to spin up Coq?

Making the bin edge computation more user-friendly is an orthogonal issue, IMO. Indeed, it would be nice when it's the case that infinite-precision arithmetic would give integer bin edges, then histogram() does as well. I have my doubts that it would solve the inconsistency problem in general.

Someone want to make a PR for this.

Unfortunately I don't have time to do this in the short term - it could be a few weeks until I have time to properly look over this. I do think it would be worth trying to fix linspace first before considering completely reverting the changes in np.histogram, given the performance implications. I think in 99.9% of real-life cases, users will not be concerned with bin edge accuracy down to 1-2 ULP, so while I think we do want to get it right, I don't think this is a critical/urgent issue, and it would make sense to take our time to get it right.

One possible compromise if we do need to revert np.histogram though - provide a new np.histogram_equal function that computed only histogram with equal bin sizes, and in that case does not return bin edges? (it would just return the y values of a histogram computed under the assumption of equal bin sizes, given a min/max and a number of bins). Would that be an acceptable solution? (I think this would cover many use cases, and this would be the high performance function).

I don't think histogram_equal() is an option. There is no histogram without the bins. Users will just use linspace() and get the same inconsistent results.

I don't think that you can change linspace() to always be consistent with your indexing method here, but if you want to try, here is a fuzz tester: https://gist.github.com/rkern/4ec8ebbb221191304a072eaed2473508

The other solution is more straightforward, and should take far less time to do properly; it just needs to be coded in Cython for performance. Roughly,

bin_edges = linspace(mn, mx, bins)
indices = np.empty(len(data), dtype=np.intp)
for i range(len(data)):
    j = (x - mn) * bins / (mx - mn)
    if j < len(bin_edges) - 1 and x >= bin_edges[j + 1]:
        j += 1
    if j > 0 and x < bin_edges[j]:
        j -= 1
    indices[i] = j

I would not mind going the cython route for later release, but first the histogram functions should be separated out into their own module. I want to do that anyway as the relevant code has expanded. What I do want to solve is the 1.11.0 regression so that I can get on with the 1.11.1 release. I sounds like a reversion would be best for that, and maybe for 1.12.0 also, as I don't want to delay the 1.12 branch unduly.

@rkern

@rkern Straight reversion is a bit of a mess due to other commits. Could you do your proposed fix in python? If the performance isn't too awful I'd like to use it. After I will separate out the histogram functions into their own file in master and set things up to use cython.

Navigation Menu

Repository navigation

np.histogram() inconsistent results #7628

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

np.histogram() inconsistent results #7628

Description

Activity

rkern commented on May 12, 2016

charris commented on May 14, 2016

rkern commented on May 15, 2016

rkern commented on May 15, 2016

rkern commented on May 15, 2016

charris commented on May 16, 2016

rkern commented on May 16, 2016

astrofrog commented on May 16, 2016

rkern commented on May 16, 2016

charris commented on May 22, 2016

astrofrog commented on May 22, 2016

rkern commented on May 22, 2016

charris commented on May 22, 2016

charris commented on May 23, 2016

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions