Fast sampling from discrete distributions

vancleve · October 15, 2018, 11:40pm

I’m writing some population simulations that make heavy use of categorical/multinomial sampling. I can see from the code that the categorical and multinomial samplers use the alias table method.

In this multinomial, the alias method is only used when n^2 < k where n is the number of trials and k is the number of possible events.

github.com

JuliaStats/Distributions.jl/blob/8bdcda68f903cf377ce37ee4540a3e4d936877e1/src/samplers/multinomial.jl#L53


      
              prob::Vector{Float64}
              alias::AliasTable
          end
          
          
MultinomialSampler(n::Int, prob::Vector{Float64}) = 
              MultinomialSampler(n, prob, AliasTable(prob))
          
          
function _rand!(s::MultinomialSampler, x::AbstractVector{T}) where T<:Real
              n = s.n
              k = length(s)
              if n^2 > k      
                  multinom_rand!(n, s.prob, x)
              else
                  # Use an alias table
                  fill!(x, zero(T))
                  a = s.alias
                  for i = 1:n
                      x[rand(a)] += 1
                  end
              end
              return x

At some point I think that @johnmyleswhite added this cutoff; I was wondering what the rationale for it was?

Also, I wonder if there are some faster methods available. For example, https://arxiv.org/pdf/1611.00532.pdf
present an algorithm (#5 in Table 1) that appears to scale better with k than the alias table method (Walker’s algorithm). I don’t know the literature well, so I thought I’d mention this in case its new to folks but am generally interested in what are the fastest methods (big surprise ).