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.

Wimmer

The binomial distribution counts the number of "successes" occurring in n trials when the trials are independent and have a fixed probability of success. No successes is one of the possibilities. Since your algorithm never generates zeroes, it is not generating a binomial distribution. – Carlcarla

@Carlcarla I meant Multinomial with one draw. Someone justly corrected the title. – Wimmer

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

Babbette

I figured it out. But thanks a bunch for teaching me about proptable! – Wimmer

Answer 2

The solution is to use sampling with weighted probabilities.

Add StatsBase package if not yet added

Pkg.add("StatsBase")

Draw samples:

using StatsBase
const probabilities = [0.3, 0.3, 0.2, 0.15, 0.05]
items = [i for i in 1:length(probabilities)]
weights = Weights(probabilities)

for i in 1:20
    println(sample(items, weights))
end

Wimmer