My mistake (though this optimisation has nothing to do with accepting the PR anyway). This lookup table method is (a little) faster for 1e6 samples, as long as you use:
```python
prob_array = np.digitize(np.linspace(0, 1, map_length), probs)
```
or
```python
prob_array = np.searchsorted(probs, np.linspace(0, 1, map_length))
```
instead of:
```python
prob_array = np.array([np.searchsorted(probs, r) for r in np.linspace(0, 1, map_length)])
```
(`searchsorted` is better as the number of classes grows, with parity around 100 classes on my system, but with your lookup table the difference is marginal.)

But for 1e8 samples the number of intermediate arrays still makes the most difference. You can remove one intermediate copy by removing `np.floor`: for non-negative data, `astype(int)` should be the same as `floor`. Because `rs` is only used here, you can remove another by using `rs *= map_length`.