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()
ISTM that this elementary functionality deserves an implementation that's as fast as it can be.