title: numpy mean of multidimensional array
body: <p>I have a multidimensional numpy array that happens to be an array of images. Why does computing the image channel mean produce different results when using the axis argument to np.mean?</p>
  
  <pre><code>&gt;&gt;&gt; X = np.array(np.random.random((9999, 128, 128, 4)) * 1e5, dtype='float32')
  &gt;&gt;&gt; X.shape
  (9999, 128, 128, 4)
  &gt;&gt;&gt; mean_by_axis = np.mean(X, axis=(0, 1, 2))
  array([ 13423.11523438,  13423.11523438,  13423.11523438,  13423.11523438], dtype=float32)
  &gt;&gt;&gt; mean = np.mean(X[:, :, :, 0])
  50001.297
  </code></pre>
  
  <p>I expect mean_by_axis[0] == mean. Why is this not the case? The same is true of the remaining axis-3 indices 1, 2 and 3. Am I misunderstanding how to use the axis argument in np.mean?</p>
  
  <p>Using numpy version '1.12.1'</p>
  
  <p>Is it possible that I am overflowing the float32 accumulator? For instance:</p>
  
  <pre><code>&gt;&gt;&gt; X = np.random.random(size=(100, 128, 128, 4))
  &gt;&gt;&gt; np.mean(X, axis=(0, 1, 2))
  array([ 0.49978557,  0.49985835,  0.50000321,  0.50015689])]
  &gt;&gt;&gt; np.mean(X[:, :, :, 0])
  0.49978556940636332
  </code></pre>
  
  <p>This looks about correct. If this is the case, why does the slice method not also overflow the accumulator and give the same result? Perhaps the slice methods uses a float64 accumulator and the axis method uses a float32 accumulator?</p>