title: Numpy average from a large masked array
body: <p>Whats the appropriate way to get the average value from a large masked array? Normally i would just call <code>.mean()</code> but for very large array's this fails for me. </p>
  
  <p>Consider creating an array of a million elements, all with value 500, like:</p>
  
  <pre><code>a = np.ones(1000000, dtype=np.int16) * 500
  </code></pre>
  
  <p>Then create a random mask and combine both in a new <code>masked array</code>:</p>
  
  <pre><code>mask = np.random.randint(0, 2, a.size)
  b = np.ma.masked_array(a, mask=mask)
  </code></pre>
  
  <p>The array <code>b</code> inherits the <code>dtype</code> and is also a <code>int16</code>.</p>
  
  <p>Getting the average value from <code>b</code> can be done in different ways, but all give the same result. But the <code>non-ma</code> functions are ignoring the mask and shouldn't be used.</p>
  
  <pre><code>print(np.average(b), np.mean(b), np.ma.average(b), np.ma.mean(b))
  (500.0, 500.0, 500.0, 500.0)
  </code></pre>
  
  <p>But if i increase the original array size from one million to ten million, the result becomes:</p>
  
  <pre><code>print(np.average(b), np.mean(b), np.ma.average(b), np.ma.mean(b))
  (500.0, -359.19365132075774, -359.19365132075774, -359.19365132075774)
  </code></pre>
  
  <p>Now only <code>np.average</code> <strong>seems</strong> correct, but as said is ignoring the mask and calculating the average over the entire array, this can be shown when changing some of the masked values with <code>b[b.mask] = 1000</code> for example. I would expect <code>np.mean</code> to do the same though.</p>
  
  <p>Casting the masked array <code>b</code> to a <code>float32</code> results in:</p>
  
  <pre><code>b = b.astype(np.float32)
  print(np.average(b), np.mean(b), np.ma.average(b), np.ma.mean(b))
  (511.18945, 510.37895680000003, 510.37895680000003, 510.37895680000003)
  </code></pre>
  
  <p>And casting the masked array <code>b</code> to a <code>float64</code> (which should be done by default according to the documentation) results in:</p>
  
  <pre><code>b = b.astype(np.float64)
  print(np.average(b), np.mean(b), np.ma.average(b), np.ma.mean(b))
  (500.0, 500.0, 500.0, 500.0)
  </code></pre>
  
  <p>So casting to a <code>float64</code> seems to work, but i would rather avoid this since it would increase the memory footprint a lot.</p>
  
  <p>And while testing some of this, i also noted that calling <code>np.ma.average</code> on the non-masked array (a) gives the correct result if the size is one million and the wrong result when its ten million, while <code>np.ma.mean</code> is correct for both sizes.</p>
  
  <p>Can someone explain the relation between <code>dtype</code> and the <code>size</code> of the array in this example? Its a bit of a mystery to me when this occurs and how to handle it properly.</p>
  
  <p>All is done in Numpy 1.8.1 on a 64bit Win 7 machine. Installed via conda.</p>
  
  <p>Here is a notebook replicating what i have done:</p>
  
  <p><a href="http://nbviewer.ipython.org/gist/RutgerK/69b60da73f464900310a" rel="nofollow">http://nbviewer.ipython.org/gist/RutgerK/69b60da73f464900310a</a></p>