title: Calculate histograms along axis
body: <p>Is there a way to calculate many histograms along an axis of an nD-array? The method I currently have uses a <code>for</code> loop to iterate over all other axes and calculate a <code>numpy.histogram()</code> for each resulting 1D array:</p>
  
  <pre><code>import numpy
  import itertools
  data = numpy.random.rand(4, 5, 6)
  
  # axis=-1, place `200001` and `[slice(None)]` on any other position to process along other axes
  out = numpy.zeros((4, 5, 200001), dtype="int64")
  indices = [
      numpy.arange(4), numpy.arange(5), [slice(None)]
  ]
  
  # Iterate over all axes, calculate histogram for each cell
  for idx in itertools.product(*indices):
      out[idx] = numpy.histogram(
          data[idx],
          bins=2 * 100000 + 1,
          range=(-100000 - 0.5, 100000 + 0.5),
      )[0]
  
  out.shape  # (4, 5, 200001)
  </code></pre>
  
  <p>Needless to say this is very slow, however I couldn't find a way to solve this using <code>numpy.histogram</code>, <code>numpy.histogram2d</code> or <code>numpy.histogramdd</code>.</p>