# Draw a random number through a probability distribution defined as an array

**URL:** <https://discourse.julialang.org/t/draw-a-random-number-through-a-probability-distribution-defined-as-an-array/129210>\
**Category:** General Usage\
**Tags:** distributions, random, sampling\
**Created:** [May 21, 2025, 9:36am UTC](https://discourse.julialang.org/t/draw-a-random-number-through-a-probability-distribution-defined-as-an-array/129210 "2025-05-21T09:36:47Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [May 21, 2025, 9:36am UTC](https://discourse.julialang.org/t/draw-a-random-number-through-a-probability-distribution-defined-as-an-array/129210/1 "2025-05-21T09:36:47Z")

</div>

If I have an NxM matrix with decimal values (normalised so that the sum is 1):

```julia
prob_dens::Matrix{Float64.Float64}

```

is there a good way to randomly draw an index `i,j` with probability corresponding to the probability distribution at that index?

---

<div class="post-metadata">

**Author:** ![Alexander-Barth](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alexander-barth/32/3692_2.png) [@Alexander-Barth](https://discourse.julialang.org/u/Alexander-Barth)\
**Post date:** [May 21, 2025, 10:08am UTC](https://discourse.julialang.org/t/draw-a-random-number-through-a-probability-distribution-defined-as-an-array/129210/2 "2025-05-21T10:08:08Z")

</div>

You can get an `i,j` candidate at random (uniform distribution) and accept it with a probability of prob\_dens[i,j].

```julia
function sample(prob_dens)
while true
    I = rand(CartesianIndices(prob_dens))
    if rand() < prob_dens[I]
        return I
    end
end
end

prob_dens = rand(4,4); prob_dens /= sum(prob_dens)
sample(prob_dens)

```

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [May 21, 2025, 10:08am UTC](https://discourse.julialang.org/t/draw-a-random-number-through-a-probability-distribution-defined-as-an-array/129210/3 "2025-05-21T10:08:51Z")

</div>

The following worked for me:

```julia
using StatsBase

sample(CartesianIndices(prob_dens), weights(vec(prob_dens)))

```

The output can be wrapped with `Tuple` to get a tuple of the indices.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [May 21, 2025, 10:29am UTC](https://discourse.julialang.org/t/draw-a-random-number-through-a-probability-distribution-defined-as-an-array/129210/4 "2025-05-21T10:29:38Z")

</div>

It’s probably best to stick to some established proper implementation, but here’s a ‘home-made’, which isn’t really tested, but shows an approach:

```julia
# Create normalized matrix
M = rand(5, 6);
M ./= sum(M)

inds = CartesianIndices(M) 
C = cumsum(vec(M)) # C is now a sort of cumulative distribution

function mysample(C, inds)
    i = searchsortedfirst(C, rand())
    return inds[i]
end

```

It uses a single `rand()` call and then looks up the corresponding element in the cumulative distribution. I’m not confident that it isn’t off-by-one or whether to use `<` or `<=` internally. But a single `rand()` is fast, and searching sorted vectors is fast.

---

<div class="post-metadata">

**Author:** ![Torkel](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/torkel/32/5030_2.png) [@Torkel](https://discourse.julialang.org/u/Torkel)\
**Post date:** [May 21, 2025, 10:37am UTC](https://discourse.julialang.org/t/draw-a-random-number-through-a-probability-distribution-defined-as-an-array/129210/5 "2025-05-21T10:37:31Z")

</div>

Thanks a lot everyone, that is very helpful!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 21, 2025, 10:40am UTC](https://discourse.julialang.org/t/draw-a-random-number-through-a-probability-distribution-defined-as-an-array/129210/6 "2025-05-21T10:40:30Z")

</div>

> [@Alexander-Barth](#):
>
> You can get an `i,j` candidate at random (uniform distribution) and accept it with a probability of prob\_dens[i,j].

This will be slow for a large array — the expected number of times it calls `rand()` is proportional to the length of the array.

> [@Dan](#):
>
> `sample(CartesianIndices(prob_dens), weights(vec(prob_dens)))`

Note that this also takes time proportional to `length(prob_dens)` (though is faster because it only generates a single random number and scans the array consecutively), because the algorithm works by [scanning a cumsum through the weights array](https://github.com/JuliaStats/StatsBase.jl/blob/dfe8945a90638de0a42d8bc7a8eb6f4328d2a358/src/sampling.jl#L586-L601) to do [inverse transform sampling](https://en.wikipedia.org/wiki/Inverse_transform_sampling).

If you need to generate many samples for a fixed weights table, it is probably faster to use AliasTables.jl, which allows you to precompute an efficient mapping from random numbers to indices via the [alias method](https://en.wikipedia.org/wiki/Alias_method):

```julia
using AliasTables
at = AliasTable(vec(prob_dens))
CartesianIndices(prob_dens)[rand(at)]

```

e.g. this is almost 200x faster for me than `sample` for `prob_dens = rand(100,200)` (not including the one-time cost to generate the `AliasTable`).