20

I want to get a list of the data contained in a histogram bin. I am using numpy, and Matplotlib. I know how to traverse the data and check the bin edges. However, I want to do this for a 2D histogram and the code to do this is rather ugly. Does numpy have any constructs to make this easier?

For the 1D case, I can use searchsorted(). But the logic is not that much better, and I don’t really want to do a binary search on each data point when I don’t have to.

Most of the nasty logic is due to the bin boundary regions. All regions have boundaries like this: [left edge, right edge). Except the last bin, which has a region like this: [left edge, right edge].

Here is some sample code for the 1D case:

import numpy as np

data = [0, 0.5, 1.5, 1.5, 1.5, 2.5, 2.5, 2.5, 3]

hist, edges = np.histogram(data, bins=3)

print 'data =', data
print 'histogram =', hist
print 'edges =', edges

getbin = 2  #0, 1, or 2

print '---'
print 'alg 1:'

#for i in range(len(data)):
for d in data:
    if d >= edges[getbin]:
        if (getbin == len(edges)-2) or d < edges[getbin+1]:
            print 'found:', d
        #end if
    #end if
#end for

print '---'
print 'alg 2:'

for d in data:
    val = np.searchsorted(edges, d, side='right')-1
    if val == getbin or val == len(edges)-1:
        print 'found:', d
    #end if
#end for

Here is some sample code for the 2D case:

import numpy as np

xdata = [0, 1.5, 1.5, 2.5, 2.5, 2.5, \
         0.5, 0.5, 0.5, 0.5, 1.5, 1.5, 1.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, \
         0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 3]
ydata = [0, 5,5, 5, 5, 5, \
         15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, 15, \
         25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 30]

xbins = 3
ybins = 3
hist2d, xedges, yedges = np.histogram2d(xdata, ydata, bins=(xbins, ybins))

print 'data2d =', zip(xdata, ydata)
print 'hist2d ='
print hist2d
print 'xedges =', xedges
print 'yedges =', yedges

getbin2d = 5  #0 through 8

print 'find data in bin #', getbin2d

xedge_i = getbin2d % xbins
yedge_i = int(getbin2d / xbins) #IMPORTANT: this is xbins

for x, y in zip(xdata, ydata):
    # x and y left edges
    if x >= xedges[xedge_i] and y >= yedges[yedge_i]:
        #x right edge
        if xedge_i == xbins-1 or x < xedges[xedge_i + 1]:
            #y right edge
            if yedge_i == ybins-1 or y < yedges[yedge_i + 1]:
                print 'found:', x, y
            #end if
        #end if
    #end if
#end for

Is there a cleaner / more efficient way to do this? It seems like numpy would have something for this.

Ben
  • 1,038
  • 1
  • 19
  • 22
  • 3
    Just out of curiosity; why do you use comments like #end if in your code? "Every pixel counts" By doing that you are ignoring the purpose of indentation. – Gökhan Sever Feb 16 '10 at 23:18
  • 3
    2 reasons. I am a C++ developer first, and a python developer second. Python's lack of braces irritates me to no end. When I have complicated code blocks with lots of varying indentation, I don't want to be counting whitespace. And I do most of my development in Emacs. By putting closing comments on code blocks, it lets me press TAB on every line and Emacs won't try to wrongly indent something. – Ben Feb 19 '10 at 02:53

2 Answers2

27

digitize, from core NumPy, will give you the index of the bin to which each value in your histogram belongs:

import numpy as NP
A = NP.random.randint(0, 10, 100)

bins = NP.array([0., 20., 40., 60., 80., 100.])

# d is an index array holding the bin id for each point in A
d = NP.digitize(A, bins)
Steven C. Howell
  • 16,902
  • 15
  • 72
  • 97
doug
  • 69,080
  • 24
  • 165
  • 199
  • This is almost perfect! If there are any numpy devs here, this function should really go in the "see also" part of the histogram docs. It is too bad that the digitize() bin logic doesn't exactly match the histogram() bin logic though. So this leads to just as awkward code as the rest of my examples above. – Ben Feb 19 '10 at 02:48
  • 3
    isn't this exactly the same as `bins.searchsorted(A, 'right')`? – Bi Rico Feb 29 '12 at 23:38
6

how about something like:

data = numpy.array([0, 0.5, 1.5, 1.5, 1.5, 2.5, 2.5, 2.5, 3])

hist, edges = numpy.histogram(data, bins=3)

for l, r in zip(edges[:-1], edges[1:]):
   print(data[(data > l) & (data < r)])

Out:

[ 0.5]
[ 1.5  1.5  1.5]
[ 2.5  2.5  2.5]

with a bit of code to handle the edge cases.

Tomerikoo
  • 18,379
  • 16
  • 47
  • 61
Autoplectic
  • 7,566
  • 30
  • 30