<p dir="auto">Calculating <code class="notranslate">std</code> on <code class="notranslate">float32</code> arrays using <code class="notranslate">axis</code> can give inaccurate results. See the example below.</p>
<h3 dir="auto">Reproducing code example:</h3>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="import numpy as np

n = 1_000_000

# Create a float32 array with two columns with values 100 and -100
# x.shape = (n, 2)
x = np.array([100.0, -100.0])[None, :].repeat(repeats=n, axis=0).astype(np.float32)

# calculating std for the column separately gives 0 as expected
print(x[:, 0].std())  # -&gt; 0.
print(x[:, 1].std())  # -&gt; 0.

# calculating std on both columns in one go using axis=0 gives
print(x.std(axis=0))  # -&gt; [1.3201232, 1.3201232]"><pre class="notranslate"><span class="pl-k">import</span> <span class="pl-s1">numpy</span> <span class="pl-k">as</span> <span class="pl-s1">np</span>

<span class="pl-s1">n</span> <span class="pl-c1">=</span> <span class="pl-c1">1_000_000</span>

<span class="pl-c"># Create a float32 array with two columns with values 100 and -100</span>
<span class="pl-c"># x.shape = (n, 2)</span>
<span class="pl-s1">x</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">array</span>([<span class="pl-c1">100.0</span>, <span class="pl-c1">-</span><span class="pl-c1">100.0</span>])[<span class="pl-c1">None</span>, :].<span class="pl-en">repeat</span>(<span class="pl-s1">repeats</span><span class="pl-c1">=</span><span class="pl-s1">n</span>, <span class="pl-s1">axis</span><span class="pl-c1">=</span><span class="pl-c1">0</span>).<span class="pl-en">astype</span>(<span class="pl-s1">np</span>.<span class="pl-s1">float32</span>)

<span class="pl-c"># calculating std for the column separately gives 0 as expected</span>
<span class="pl-en">print</span>(<span class="pl-s1">x</span>[:, <span class="pl-c1">0</span>].<span class="pl-en">std</span>())  <span class="pl-c"># -&gt; 0.</span>
<span class="pl-en">print</span>(<span class="pl-s1">x</span>[:, <span class="pl-c1">1</span>].<span class="pl-en">std</span>())  <span class="pl-c"># -&gt; 0.</span>

<span class="pl-c"># calculating std on both columns in one go using axis=0 gives</span>
<span class="pl-en">print</span>(<span class="pl-s1">x</span>.<span class="pl-en">std</span>(<span class="pl-s1">axis</span><span class="pl-c1">=</span><span class="pl-c1">0</span>))  <span class="pl-c"># -&gt; [1.3201232, 1.3201232]</span></pre></div>
<p dir="auto">The magnitude of the error becomes larger with larger <code class="notranslate">n</code> and larger values for the columns.<br>
It seems like this loss of accuracy should be avoidable.</p>
<h3 dir="auto">Numpy/Python version information:</h3>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; import sys, numpy; print(numpy.__version__, sys.version)
1.17.2 3.7.4 (default, Sep  7 2019, 18:27:02)
[Clang 10.0.1 (clang-1001.0.46.4)]"><pre class="notranslate"><code class="notranslate">&gt;&gt;&gt; import sys, numpy; print(numpy.__version__, sys.version)
1.17.2 3.7.4 (default, Sep  7 2019, 18:27:02)
[Clang 10.0.1 (clang-1001.0.46.4)]
</code></pre></div>
<p dir="auto">A mean over an array containing only ones should obviously return only ones. However,</p>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; import numpy as np
&gt;&gt;&gt; test_array = np.ones((10000000, 4, 15), dtype = np.float32)
&gt;&gt;&gt; print(test_array.mean(axis=(0,1)))
[ 0.4194304  0.4194304  0.4194304  0.4194304  0.4194304  0.4194304
  0.4194304  0.4194304  0.4194304  0.4194304  0.4194304  0.4194304
  0.4194304  0.4194304  0.4194304]"><pre class="notranslate"><span class="pl-c1">&gt;&gt;</span><span class="pl-c1">&gt;</span> <span class="pl-k">import</span> <span class="pl-s1">numpy</span> <span class="pl-k">as</span> <span class="pl-s1">np</span>
<span class="pl-c1">&gt;&gt;</span><span class="pl-c1">&gt;</span> <span class="pl-s1">test_array</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-en">ones</span>((<span class="pl-c1">10000000</span>, <span class="pl-c1">4</span>, <span class="pl-c1">15</span>), <span class="pl-s1">dtype</span> <span class="pl-c1">=</span> <span class="pl-s1">np</span>.<span class="pl-s1">float32</span>)
<span class="pl-c1">&gt;&gt;</span><span class="pl-c1">&gt;</span> <span class="pl-en">print</span>(<span class="pl-s1">test_array</span>.<span class="pl-en">mean</span>(<span class="pl-s1">axis</span><span class="pl-c1">=</span>(<span class="pl-c1">0</span>,<span class="pl-c1">1</span>)))
[ <span class="pl-c1">0.4194304</span>  <span class="pl-c1">0.4194304</span>  <span class="pl-c1">0.4194304</span>  <span class="pl-c1">0.4194304</span>  <span class="pl-c1">0.4194304</span>  <span class="pl-c1">0.4194304</span>
  <span class="pl-c1">0.4194304</span>  <span class="pl-c1">0.4194304</span>  <span class="pl-c1">0.4194304</span>  <span class="pl-c1">0.4194304</span>  <span class="pl-c1">0.4194304</span>  <span class="pl-c1">0.4194304</span>
  <span class="pl-c1">0.4194304</span>  <span class="pl-c1">0.4194304</span>  <span class="pl-c1">0.4194304</span>]</pre></div>
<p dir="auto">This returns the correct result:</p>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; print(test_array.mean(axis=0).mean(axis=0))
[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]"><pre class="notranslate"><span class="pl-c1">&gt;&gt;</span><span class="pl-c1">&gt;</span> <span class="pl-en">print</span>(<span class="pl-s1">test_array</span>.<span class="pl-en">mean</span>(<span class="pl-s1">axis</span><span class="pl-c1">=</span><span class="pl-c1">0</span>).<span class="pl-en">mean</span>(<span class="pl-s1">axis</span><span class="pl-c1">=</span><span class="pl-c1">0</span>))
[ <span class="pl-c1">1.</span>  <span class="pl-c1">1.</span>  <span class="pl-c1">1.</span>  <span class="pl-c1">1.</span>  <span class="pl-c1">1.</span>  <span class="pl-c1">1.</span>  <span class="pl-c1">1.</span>  <span class="pl-c1">1.</span>  <span class="pl-c1">1.</span>  <span class="pl-c1">1.</span>  <span class="pl-c1">1.</span>  <span class="pl-c1">1.</span>  <span class="pl-c1">1.</span>  <span class="pl-c1">1.</span>  <span class="pl-c1">1.</span>]</pre></div>
<p dir="auto">I guess the reason for this problem is some overflow, because it does not appear when I use a test array with dtype = np.float64. I would expect numpy to either give the correct result or to at least give a warning whenever such an overflow happens.</p>
<p dir="auto">surprisingly mean along all axis gives the correct result again:</p>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="&gt;&gt;&gt; print(test_array.mean(axis=(0,1,2)))
1.0"><pre class="notranslate"><span class="pl-c1">&gt;&gt;</span><span class="pl-c1">&gt;</span> <span class="pl-en">print</span>(<span class="pl-s1">test_array</span>.<span class="pl-en">mean</span>(<span class="pl-s1">axis</span><span class="pl-c1">=</span>(<span class="pl-c1">0</span>,<span class="pl-c1">1</span>,<span class="pl-c1">2</span>)))
<span class="pl-c1">1.0</span></pre></div>
<p dir="auto">(I used numpy-1.12.1 with python 3.5 on Ubuntu 16)</p>