For n unequal weights and k selections, sample selection with the inverse-cdf method is O(k log₂ n).  Using the alias method, it improves to O(k).  The proportionally constants also favor the alias method so that if the set up times were the same, the alias method would always win (even when n=2).

However, the set up times are not the same.  For the inverse-cdf method, set up is O(1) if cum_weights are given; otherwise, it is O(n) with a fast loop.  The setup time for the alias method is also O(n) but is proportionally much slower.

So, there would need to be a method selection heuristic based on the best trade-off between setup time and sample selection time.

Both methods make k calls to random().

See: https://en.wikipedia.org/wiki/Alias_method

Notes on the attached draft implementation:

* Needs to add back the error checking code.

* Need a better method selection heuristic.

* The alias table K defaults to the original index
  so that there is always a valid selection even
  if there are small rounding errors.

* The condition for the aliasing loop is designed
  to have an early-out when the remaining blocks
  all have equal weights.  Also, the loop condition
  makes sure that the pops never fail even if there
  are small rounding errors when partitioning
  oversized bins or if the sum of weights isn't
  exactly 1.0.

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.