Well, numpy is not overly fancy about how to calculate means, there is a thing that you get better then naive precision in some cases (which kicks in for the full array here), see also for example gh-8116

Ok. So I get why a naive summation with float32 would that, but then it was quite confusing to see it succeed without axis selection.
As I understand from your comment this is due to pairwise summation in some of the cases (?)
Is this strange behavior what one should expect?
Why not implement pairwise summation for all situations?
And (maybe this naive) - is there a way to warn about such drastic round errors without a huge performance hit?
(these type of bugs are usually hard to pin down)

The reason is memory layout and thus speed. Doing the (mostly) pairwise summation in numpy with the typical machinery only works reasonably along a single axis (where it comes with no performance loss at all). But only when summing the fast axis it is feasable to do this, because otherwise others would be complaining about massive performance drops.

I agree that there should be more documentation on this, heck I was even hesitant when we first put this in.... It would also be nice to have more stable summations in general....

I'm glad to see #9393 that will add something in the documentation. But this can be a serious issue that can go unnoticed. A warning would be even more welcome. A fix is of course the idea situation.

My recommended workaround for people running into this issue is to add ,dtype=np.float64 to the np.mean() and np.std() calls. Doing that in the above code yields the following result:

[127.48983901 127.50801956 127.49455946]
[127.48983901 127.50801956 127.49455946]
127.4976

Another minimalist example that illustrates this issue:

a=np.random.rand(30*1000*1000,2).astype(np.float32)*.01+3.9 # 30 million pairs
print('Expected:')
print(' mean: ', np.mean(a,axis=0,dtype=np.float64), '  (analytical = 3.905)')
print('  std: ', np.std(a,axis=0,dtype=np.float64), '  (analytical = .01/sqrt(12) = 0.00288...)')
print('Instead, you get:')
print(' mean: ', np.mean(a,axis=0))
print('  std: ', np.std(a,axis=0))

Outputs:

Expected:
 mean:  [3.9050007  3.90499964]   (analytical = 3.905)
  std:  [0.00288722 0.00288692]   (analytical = .01/sqrt(12) = 0.00288...)
Instead, you get:
 mean:  [2.236962 2.236962]
  std:  [1.4956477 1.4956477]

I expect this situation to be encountered a lot and without noticing in neural network data normalization code where statistics over large training data gets computed. People tend to go for single (or half) rather than double precision because of performance and memory savings, and the fact that successful neural network training rarely requires double precision.

I suppose the warning could be made length dependent.