title: Numpy histogram on multi-dimensional array
body: <p>given an np.array of shape <code>(n_days, n_lat, n_lon)</code>, I'd like to compute a histogram with fixed bins for each lat-lon cell (ie the distribution of daily values). </p>
  
  <p>A simple solution to the problem is to loop over the cells and invoke <code>np.histogram</code> for each cell::</p>
  
  <pre><code>bins = np.linspace(0, 1.0, 10)
  B = np.rand(n_days, n_lat, n_lon)
  H = np.zeros((n_bins, n_lat, n_lon), dtype=np.int32)
  for lat in range(n_lat):
      for lon in range(n_lon):
          H[:, lat, lon] = np.histogram(A[:, lat, lon], bins=bins)[0]
  # note: code not tested
  </code></pre>
  
  <p>but this is quite slow. Is there a more efficient solution that does not involve a loop?</p>
  
  <p>I looked into <code>np.searchsorted</code> to get the bin indices for each value in <code>B</code> and then use fancy indexing to update <code>H</code>::</p>
  
  <pre><code>bin_indices = bins.searchsorted(B)
  H[bin_indices.ravel(), idx[0], idx[1]] += 1  # where idx is a index grid given by np.indices
  # note: code not tested
  </code></pre>
  
  <p>but this does not work because the in-place add operator (+=) doesn't seem to support multiple updates of the same cell.</p>
  
  <p>thx, 
   Peter</p>