but this is bound to be inefficient as soon as the vector of probabilities gets large, especially if you want to draw multiple samples.

Have I overlooked something or should this be added?

I think you misunderstand the point of multinomial distributions.

I'm afraid the multiple samples were an afterthought and a false lead. My main use case would be consecutive samples of _different_ categorical distributions, for which abusing multinomials seems wasteful (as they have to allocate a large vector each time). Similarly for the alternative spelling

import numpy, random a = numpy.array([.5, .3, .2]) (a.cumsum()-random.random() >= 0).nonzero()[0][0]

1

ISTM that this elementary functionality deserves an implementation that's as fast as it can be.

- Hagen