DOC: Document higher precision in mean/sum/add.reduce along fast axis due to pairwise sum. #9393
Description
As mentioned in this stackoverflow question, calling np.mean with the axis argument seems to introduce significantly lower numerical stability. In particular,
>>> import numpy as np
>>> X = np.random.rand(9999, 128, 128, 4).astype('float32')
>>> X.shape
>>> np.mean(X, axis=(0, 1, 2))
array([ 0.10241024, 0.10241024, 0.10241024, 0.10241024], dtype=float32) # should be 0.5
>>> np.mean(X[:, :, :, 0])
0.50000387whereas calling with full 64 bit precision gives correct values
>>> np.mean(X.astype('float64'), axis=(0, 1, 2))
array([ 0.50000323, 0.50004907, 0.50003198, 0.49999848])
>>> np.mean(X[:, :, :, 0].astype('float64'))
0.50000323305421812Since the values are correct without the axis argument, the should likely be classified as a bug.
Edit:
This is on windows 7 with numpy 1.12.0, but it seems others are experiencing this issue as well.