Distributions.jl using a custom random number generator

bernhard · July 28, 2017, 9:39am

Is there any way to use a distribution from Distributions.jl with a custom RNG?

I know how to sample in a reproducible manner from a uniform distribution (see x1 and x2 in the code below).
But Distrubutions.jl seems to rely on Base.Random.GLOBAL_RNG which I do not want to use

I read this post, but was unable to solve my issue

I found a Julia implementation for Binomial (in LightGraphs), which I might just use, but I would prefer to have this working for any distribution.

github.com

sbromberger/LightGraphs.jl/blob/1d4c341fdad3cfb6954361c107c7a5ca32903ae2/src/generators/randgraphs.jl#L41


      
                  dest = rand(rng, one(T):tnv)
                  source != dest && add_edge!(g, source, dest)
              end
              return g
          end
          
          DiGraph(nv::T, ne::Integer; seed::Int = -1) where T<:Integer =
              DiGraph{Int}(nv, ne, seed=seed)
          
          """
              randbn(n, p, seed=-1)
          
          Return a binomally-distribted random number with parameters `n` and `p` and optional `seed`.
          
          ### References
          - "Non-Uniform Random Variate Generation," Luc Devroye, p. 522. Retrieved via http://www.eirene.de/Devroye.pdf.
          - http://stackoverflow.com/questions/23561551/a-efficient-binomial-random-number-generator-code-in-java
          """
          function randbn(n::Integer, p::Real, seed::Integer=-1)
              rng = getRNG(seed)
              log_q = log(1.0 - p)

using RandomNumbers, Distributions

const RNG_XOR=Xorshifts.Xorshift1024Plus(0x1234567890abcdef)
const RNG_MERSENNE=MersenneTwister(2001)
const RNG_SELECTED=RNG_MERSENNE

x1=rand(RNG_XOR) #this works fine, x1 is reproducible each time I initialize my seed & custom RNG, I get the same value for x1
x2=rand(RNG_MERSENNE) #this works fine as well

# this line throws an error (I was hoping this would solve my problem...)
# Base.Random.GLOBAL_RNG=RNG_SELECTED

some_distr=Binomial(200,0.1)
rand(some_distr) #this 'value' is only reproducible if i run srand(some_seed) beforehand, but I want to use a different RNG

evanfields · July 28, 2017, 7:56pm

As I understand it, this isn’t currently possible for all distributions. Distributions.jl uses Rmath for various complicated functions, and so Julia RNG objects can’t be passed to control random number generation.

If this isn’t the case I’d love to hear it. Sampling from Distributions.jl with a specified RNG is a much-desired feature for me, and I think the reliance on Rmath is a (totally understandable) wart in the Julia ecosystem.

bernhard · July 29, 2017, 6:12am

Thank you. I saw that RMath was called.
Luckily Binomial, Bernoulli and Normal suffice for my use case (for now). Thus I can remain within Julia and rely on a custome RNG for those.