title: Calculate column-wise mean of rows of a numpy matrix selected using values from another array
body: <p>Given a numpy array <code>M</code> of shape <code>(r, c)</code> e.g.</p>
  <pre><code>M = np.array([[1, 2, 3],
                [4, 5, 6],
                [7, 8, 9],
                [10, 11, 12],
                [13, 14, 15]])  # r = 5; c = 3
  </code></pre>
  <p>and a one-dimensional array <code>a</code> of length <code>r</code>, containing integers that can vary from 0 to <code>k-1</code>, e.g.</p>
  <pre><code>a = np.array([0, 0, 2, 1, 0])  # k = 4
  </code></pre>
  <p>I want to use the values in <code>a</code> to select rows from <code>M</code> to get an intermediate result like this:</p>
  <pre><code>array([
         [[1, 2, 3], [4, 5, 6], [13, 14, 15]]  # rows of M where a == 0
         [[10, 11, 12]],                       # rows of M where a == 1
         [[7, 8, 9]]                           # rows of M where a == 2
         []                                    # rows of M where a == 3 (there are none)
        ]) 
  </code></pre>
  <p>(I don't need this intermediate array, but am just showing it for illustration.) The returned result would be a <code>(k, c)</code> array with the column-wise means from this array:</p>
  <pre><code>array([[ 6.,  7.,  8.],   # means where a == 0
     [10., 11., 12.],       # means where a == 1
     [ 7.,  8.,  9.],       # etc.
     [nan, nan, nan]])
  </code></pre>
  <p>I can do this with</p>
  <pre><code>np.array([M[a == i].mean(axis=0) for i in range(k)])
  </code></pre>
  <p>but is there a way (hopefully faster for large <code>r</code> and <code>k</code>) that purely uses numpy methods instead of using a for loop to make a list (which will then have to be converted back to an array)?</p>