I also looked at another method using binomial variates but couldn't get it to run faster than the alias method:

    def choices(population, weights, *, k=1):
        r = 1.0
        n = k
        selections = []
        for elem, p in zip(population, weights):
            v = binomial_variate(n, p / r)
            selections += [elem] * v
            n -= v
            r -= p
        shuffle(selections)
        return selections

The shuffle step took as much time as the alias method.  Also, the binomial variate was hard to compute quickly and without overflow/underflow issues for large inputs.