Sampling from a discrete distribution in Julia

Question

I would like to repeatedly sample from a discrete distribution so as to get a single number.

The following is some code implementing what I am looking for:

const probabilities = [0.3, 0.3, 0.2, 0.15, 0.05]
const cummulative_probabilities = cumsum(probabilities)

function pickone(cummulative_probabilities)
    n = length(cummulative_probabilities)
    i = 1
    r = rand()
    while r >= cummulative_probabilities[i] && i<n 
        i+=1
    end
    return i
end 

for i in 1:20
    println(pickone(cummulative_probabilities))
end

The proposed alternative of using Distributions does not cut it, as the closest I can get is the following code:

using Random
using Distributions

const probabilities = [0.3, 0.3, 0.2, 0.15, 0.05]
mnd = Multinomial(1, probabilities)

for i in 1:20
    println(rand(mnd))
end

Alas, in this case, rand returns an entire vector with a single 1, the rest being zeros.

Answer 1

You seem to be looking for weighted sampling, which is implemented by sample in StatsBase.jl:

julia> using StatsBase

julia> using FreqTables

julia> 

julia> proptable([pickone(cummulative_probabilities) for _ in 1:10^7])
5-element Named Array{Float64,1}
Dim1  │ 
──────┼──────────
1     │  0.300094
2     │       0.3
3     │  0.199871
4     │  0.150075
5     │ 0.0499595

julia> proptable([sample(Weights(probabilities)) for _ in 1:10^7])
5-element Named Array{Float64,1}
Dim1  │ 
──────┼──────────
1     │   0.29987
2     │  0.300086
3     │  0.199956
4     │  0.150184
5     │ 0.0499035