Saving Random Number State for StatsFuns' distributions Gamma Beta

New to Julia
1

On Stackoverflow there is a post showing how to save the current state of the base random number generator, so that one can go back and produce the same set of random numbers from a set point that is already so many draws in (ie it is not a simple set the seed response).

function reset_global_rng(rng_state)
Base.Random.GLOBAL_RNG.seed = rng_state.seed
Base.Random.GLOBAL_RNG.state = rng_state.state
Base.Random.GLOBAL_RNG.vals = rng_state.vals
Base.Random.GLOBAL_RNG.idx = rng_state.idx
end

using Distributions
rs = deepcopy(Base.Random.GLOBAL_RNG)
println(rand(5))
reset_global_rng(rs)
println(rand(5))

This does not work for distributions such as Beta and Gamma which utilise the StatsFuns package random number generator instead of the Base one.
I’ve tried to write a function similar to the one above but without success. Can anyone please help (so r1 is identical to r2) in code below - the deepcopy fails so almost certain the function will fail too).

function reset_statsfuns_rng(rng_state)
StatsFuns.Random.GLOBAL_RNG.seed = rng_state.seed
StatsFuns.Random.GLOBAL_RNG.state = rng_state.state
StatsFuns.Random.GLOBAL_RNG.vals = rng_state.vals
StatsFuns.Random.GLOBAL_RNG.idx = rng_state.idx
end

rs = deepcopy(StatsFuns.Random.GLOBAL_RNG)
r1 = rand(Gamma(3,11),5)
reset_statsfuns_rng(rs)
r2 = rand(Gamma(3,11),5)

Thanks.

2

  1. Please format your code using backticks (```).

  2. Are you using this library?

AFAICT it does not use random numbers.

  1. I would recommend using

which AFAIK uses the global rng.

3

Thanks for the backticks tip. The Distributions package calls the StatsFuns package for generating certain distributions (not uniform or normal). You can see this from the code source under the distributions for Gamma and Beta on the Distributions docs. So my question (unfortunately) still stands.

4

AFAICT Distributions.jl uses the global RNG.

Also, please produce a minimal working example that at least runs and demonstrates your problem. The above code does not run, since there is no such thing as StatsFuns.Random.GLOBAL_RNG.

5

GitHub allows you to capture line numbers in source code as links so tu ou can post a direct link to the line you are talking about. Otherwise it’s not clear what you mean.

6

As an aside note, for MersenneTwister it would be cleaner/more future-proof to use copy!/copy rather than copying the fields manually; the internals will change eventually. E.g. rng_state = copy(Base.Random.GLOBAL_RNG); ...; copy!(Base.Random.GLOBAL_RNG, rng_state).

1 Like
7

Code below hopefully demonstrates what I want to do. It works perfectly as is for the normal distribution. It does not work correctly for the gamma distribution (as shown in very bottom code) - reasons given in first message. If could someone could provide actual solution code I’d be very grateful.

using Distributions

function reset_global_rng(rng_state)
    Base.Random.GLOBAL_RNG.seed = rng_state.seed
    Base.Random.GLOBAL_RNG.state = rng_state.state
    Base.Random.GLOBAL_RNG.vals = rng_state.vals
    Base.Random.GLOBAL_RNG.idx = rng_state.idx
end

srand(1234) # set seed
# Now generate 100,000 Normal rands
r = rand(Normal(3,11),10000) 
#Now from herein my model can take two different paths:store generator
rs = deepcopy(Base.Random.GLOBAL_RNG)
# Now take  path 1
path1 = rand(Normal(5,10),2)
# path2 has higher mean but more uncertainty but want to minimise
#monte carlo sampling error, hence desire to reset generator
reset_global_rng(rs)
path2 = rand(Normal(6,12),2)
#To prove the resetting worked do path3 below with same parameters as path1
reset_global_rng(rs)
path3 = rand(Normal(5,10),2)
println(path1==path3) # true

Below does not work correctly for gamma (or beta)

srand(1234) # set seed
# Now generate 100,000 Normal rands
r = rand(Normal(3,11),10000) 
#Now from herein my model can take two different paths and I want to traverse both
# store generator value
rs = deepcopy(Base.Random.GLOBAL_RNG)
# Now take  path 1
path1 = rand(Gamma(5,10),2)
# path2 has higher mean but more uncertainty but want to minimise
#monte carlo sampling error, hence desire to reset generator
reset_global_rng(rs)
path2 = rand(Gamma(6,12),2)
#To prove the resetting hasn't worked:
reset_global_rng(rs)
path3 = rand(Gamma(5,10),2)
println(path1==path3) # false
8

This works for me:

using Distributions

random_seed = UInt32[0x189bde00, 0xebb1fff3, 0x8ea4d5a7, 0x856e92e5]

srand(random_seed)

path1 = rand(Gamma(5,10),2)
path2 = rand(Gamma(6,12),2)

srand(random_seed)

path3 = rand(Gamma(5,10),2)
println(path1 == path3 ? "equal" : "unequal") # prints equal
9

… but so does your original example. Are you sure you are running this exactly as above? Please post a gist of a single julia source code file, that you run with

julia path/to/file.jl

and version information (for Julia, and at least Pkg.status("Distributions")).

10

Thanks for this - I don’t want to manually set a seed value and then later call the same value back (because I want to run some common random sequences after the seed is set). But from your code it is clear that if I can extract the current four element "seed "value of the pertinent random number at the time I want to pause (where my model takes two paths) then my problems are solved - simply reapply it using your code. Obviously though “current state” is better than “seed” as a description of this four vector “seed”. Any thoughts? Many Thanks

11

This was extracted from some RNG state as Base.Random.GLOBAL_RNG.seed.

I suspect there is an orthogonal issue here: your example works for me, while it apparently did not work for you. The steps I suggested would help investigate why.

12

It’s still not clear for me why you don’t want to use srand: nothing prevents you from using srand(1234) at the start of your code snippet, and then to do something like srand(5678) before each path?

13

Sorry, are you saying that the very last line on my last posted code returned True (ie my Gamma example)? It returns False on mine (which is my whole problem). If that is the case then it is very strange

14

I don’t want to do this because the example I have given is very simple and my true model wants to do such steps at many different points, often after each other, in unpredictable ways. I do not want to keep going back and having to run past iterations that are not needed if I can keep on storing states.

15

Strange indeed, that’s why I asked for further information. Which you are not providing, so I cannot help you any further.

16

I’m using JulioPro 6.0.01 (the Juno setup using Atom - installed about a month ago). On Windows 7. Distributions package is 0.14.2. The code sample I gave was a straight cut and paste from within Juno.

Can anyone else please comment on what results they get from running my code - is the very last expression from my second code posting True or False?