Implementing a custom Sampler in Distributions.jl

Specific Domains Statistics
1

I am trying to implement my own Sampler for a scaled and shifted Beta distribution with support in some arbitrary range [a, b]:

using Distributions: Beta, Sampleable, Continuous, Univariate

struct ScaledShiftedBetaSampler <: Sampleable{Univariate, Continuous}
    distribution::Beta
    a::Float64
    b::Float64
end

function Base.rand(d::ScaledShiftedBetaSampler)
    sample = rand(d.distribution)
    return sample * (d.b - d.a) + d.a
end

According to the Distributions.jl docs it is sufficient to implement the function rand(d::ScaledShiftedBetaSampler) that returns a single random number while vectorised versions are already predefined:

The package already implements a vectorized version of rand! and rand that repeatedly calls the he scalar version to generate multiple samples.

Drawing a single sample works fine:

julia> using Random: seed!; seed!(42)
Random.MersenneTwister(UInt32[0x0000002a], ...

julia> splr = ScaledShiftedBetaSampler(Beta(2,2), 3.0, 6.0)
ScaledShiftedBetaSampler(Beta{Float64}(α=2.0, β=2.0), 3.0, 6.0)

julia> rand(splr)
4.106660698978913

julia> using Plots: histogram; histogram([rand(splr) for _ in 1:100000])

grafik
grafik766×515 13.5 KB

However, drawing several samples throws an error:

julia> rand(splr, 100000)
ERROR: ArgumentError: Sampler for this object is not defined
Stacktrace:
 [1] Random.Sampler(::Type{Random.MersenneTwister}, ::Random.SamplerTrivial{ScaledShiftedBetaSampler,Float64}, ::Val{1}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Random/src/Random.jl:145
 [2] Random.Sampler(::Random._GLOBAL_RNG, ::Random.SamplerTrivial{ScaledShiftedBetaSampler,Float64}, ::Val{1}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Random/src/Random.jl:139
 [3] rand(::Random._GLOBAL_RNG, ::Random.SamplerTrivial{ScaledShiftedBetaSampler,Float64}) at /Users/julia/buildbot/worker/package_macos64/build/usr/share/julia/stdlib/v1.4/Random/src/Random.jl:253 (repeats 2 times)
 [4] rand! at /Users/Nils/.julia/packages/Distributions/RAeyY/src/univariates.jl:165 [inlined]
 [5] rand at /Users/Nils/.julia/packages/Distributions/RAeyY/src/univariates.jl:158 [inlined]
 [6] rand(::ScaledShiftedBetaSampler, ::Int64) at /Users/Nils/.julia/packages/Distributions/RAeyY/src/genericrand.jl:24
 [7] top-level scope at REPL[12]:1

From my understanding of the docs, I was expecting this to work. What am I missing?

My versions:

julia> import Pkg; Pkg.status("Distributions")
Status `~/.julia/environments/v1.4/Project.toml`
  [31c24e10] Distributions v0.23.4

julia> versioninfo()
Julia Version 1.4.2
Commit 44fa15b150* (2020-05-23 18:35 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin18.7.0)
  CPU: Intel(R) Core(TM) i5-5250U CPU @ 1.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-8.0.1 (ORCJIT, broadwell)
1 Like
2

I think the documentation is outdated. Most of the distributions in the package implement rand(rng::AbstractRNG, s::T) for their type T.

Try implementing only this:

function Base.rand(rng::AbstractRNG, d::ScaledShiftedBetaSampler)
    sample = rand(rng, d.distribution)
    return sample * (d.b - d.a) + d.a
end

I think the rest of the machinery you want will delegate to this method, including methods that use the global rng.

3

Just as a heads up, Distributions already provides this functionality, if a bit covertly.

3 Likes
4

@contradict Yes, it seems the docs are not up to date. Defining Base.rand(rng::AbstractRNG, d::ScaledShiftedBetaSampler) does the trick, thank you very much!

@johnczito LocationScale does what I want and seems like the best option here, thanks for pointing me to it!

Problem solved :+1:

5

Not quite, until the documentation has been fixed :wink: Would you mind filing an issue with Distributions.jl?

1 Like
6

You are right — I created an issue.

3 Likes