title: Fast random weighted selection across all rows of a stochastic matrix
body: <p><code>numpy.random.choice</code> allows for weighted selection from a vector, i.e.</p>
  
  <pre class="lang-py prettyprint-override"><code>arr = numpy.array([1, 2, 3])
  weights = numpy.array([0.2, 0.5, 0.3])
  choice = numpy.random.choice(arr, p=weights) 
  </code></pre>
  
  <p>selects 1 with probability 0.2, 2 with probability 0.5, and 3 with probability 0.3.</p>
  
  <p>What if we wanted to do this quickly in a vectorized fashion for a 2D array (matrix) for which each of the rows are a vector of probabilities? That is, we want a vector of choices from a stochastic matrix? This is the super slow way:</p>
  
  <pre class="lang-py prettyprint-override"><code>import numpy as np
  
  m = 10
  n = 100 # Or some very large number
  
  items = np.arange(m)
  prob_weights = np.random.rand(m, n)
  prob_matrix = prob_weights / prob_weights.sum(axis=0, keepdims=True)
  
  choices = np.zeros((n,))
  # This is slow, because of the loop in Python
  for i in range(n):
      choices[i] = np.random.choice(items, p=prob_matrix[:,i])
  </code></pre>
  
  <p><code>print(choices)</code>:</p>
  
  <pre><code>array([ 4.,  7.,  8.,  1.,  0.,  4.,  3.,  7.,  1.,  5.,  7.,  5.,  3.,
          1.,  9.,  1.,  1.,  5.,  9.,  8.,  2.,  3.,  2.,  6.,  4.,  3.,
          8.,  4.,  1.,  1.,  4.,  0.,  1.,  8.,  5.,  3.,  9.,  9.,  6.,
          5.,  4.,  8.,  4.,  2.,  4.,  0.,  3.,  1.,  2.,  5.,  9.,  3.,
          9.,  9.,  7.,  9.,  3.,  9.,  4.,  8.,  8.,  7.,  6.,  4.,  6.,
          7.,  9.,  5.,  0.,  6.,  1.,  3.,  3.,  2.,  4.,  7.,  0.,  6.,
          3.,  5.,  8.,  0.,  8.,  3.,  4.,  5.,  2.,  2.,  1.,  1.,  9.,
          9.,  4.,  3.,  3.,  2.,  8.,  0.,  6.,  1.])
  </code></pre>
  
  <p><a href="https://stackoverflow.com/q/24140114/586086">This post</a> suggests that <code>cumsum</code> and <code>bisect</code> could be a potential approach, and is fast. But while <code>numpy.cumsum(arr, axis=1)</code> can do this along one axis of a numpy array, the <a href="https://docs.python.org/2/library/bisect.html" rel="noreferrer"><code>bisect.bisect</code></a> function only works on a single array at a time. Similarly, <a href="http://docs.scipy.org/doc/numpy/reference/generated/numpy.searchsorted.html" rel="noreferrer"><code>numpy.searchsorted</code></a> only works on 1D arrays as well.</p>
  
  <p>Is there a quick way to do this using only vectorized operations?</p>