Hi All:
In the following code, we are using StableRNGs. However, the results are not reproducing. Can anybody help us? So that we can reproduce the results.
using Turing, StableRNGs
@model function example(x)
m ~ Normal(0, 1)
x ~ Normal(m, sqrt(1))
end
res = sample(StableRNG(123), example(1), NUTS(), 1000)
In what sense? What are you changing between runs?
I’m not really sure what to look for, but I tried running that code twice and the printed output looked the same:
julia> res = sample(StableRNG(123), example(1), NUTS(), 1000)
┌ Info: Found initial step size
└ ϵ = 1.6
Sampling 100%|███████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:03
Chains MCMC chain (1000×13×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 4.44 seconds
Compute duration = 4.44 seconds
parameters = m
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ess_per_sec
Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64
m 0.4905 0.6835 0.0216 0.0314 490.8002 1.0054 110.4908
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
m -0.8832 0.0403 0.4865 0.9456 1.7660
julia> res2 = sample(StableRNG(123), example(1), NUTS(), 1000)
┌ Info: Found initial step size
└ ϵ = 1.6
Sampling 100%|███████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:00
Chains MCMC chain (1000×13×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 0.07 seconds
Compute duration = 0.07 seconds
parameters = m
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ess_per_sec
Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64
m 0.4905 0.6835 0.0216 0.0314 490.8002 1.0054 6912.6784
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
m -0.8832 0.0403 0.4865 0.9456 1.7660
(jl_XIOShu) pkg> st StableRNGs Turing
Status `/private/var/folders/gj/l9rbktlj6qndlnz1nk164d280000gn/T/jl_XIOShu/Project.toml`
[860ef19b] StableRNGs v1.0.0
[fce5fe82] Turing v0.23.1
julia> versioninfo()
Julia Version 1.8.3
Commit 0434deb161e (2022-11-14 20:14 UTC)
Platform Info:
OS: macOS (arm64-apple-darwin21.3.0)
CPU: 8 × Apple M1
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-13.0.1 (ORCJIT, apple-m1)
Threads: 4 on 4 virtual cores
Environment:
JULIA_NUM_THREADS = 4
JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
I would move to single-threaded as a first step to isolate non-reproducibility.
Thank you for showing interest in my queries. This is exactly what I did. Please see below.
julia> using Turing, StableRNGs
julia> @model function example(x)
m ~ Normal(0, 1)
x ~ Normal(m, sqrt(1))
end
example (generic function with 2 methods)
julia> res1 = sample(StableRNG(123), example(1), NUTS(), 1000)
┌ Info: Found initial step size
└ ϵ = 3.2
Sampling 100%|██████████████████████████████████████████| Time: 0:00:03
Chains MCMC chain (1000×13×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 4.36 seconds
Compute duration = 4.36 seconds
parameters = m
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat e ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
m 0.5167 0.7029 0.0222 0.0301 533.6340 0.9999 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
m -0.8081 0.0530 0.5059 0.9911 1.8123
julia> res2 = sample(StableRNG(123), example(1), NUTS(), 1000)
┌ Info: Found initial step size
└ ϵ = 1.6
Sampling 100%|██████████████████████████████████████████| Time: 0:00:00
Chains MCMC chain (1000×13×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 0.04 seconds
Compute duration = 0.04 seconds
parameters = m
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat e ⋯
Symbol Float64 Float64 Float64 Float64 Float64 Float64 ⋯
m 0.4905 0.6835 0.0216 0.0314 490.8002 1.0054 ⋯
1 column omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
m -0.8832 0.0403 0.4865 0.9456 1.7660
Clearly, two different output - for two different run.
Thank you so much for showing interest in my queries. Here is the version of my system. Do you think updating Julia will help me?
julia> versioninfo()
Julia Version 1.8.2
Commit 36034abf260 (2022-09-29 15:21 UTC)
Platform Info:
OS: macOS (arm64-apple-darwin21.3.0)
CPU: 8 × Apple M1
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-13.0.1 (ORCJIT, apple-m1)
Threads: 1 on 8 virtual cores
Thanks Dan for your reply. Can you please tell me how can I ensure that I am running it on single-thread. Here is the version of my systems
julia> versioninfo()
Julia Version 1.8.2
Commit 36034abf260 (2022-09-29 15:21 UTC)
Platform Info:
OS: macOS (arm64-apple-darwin21.3.0)
CPU: 8 × Apple M1
WORD_SIZE: 64
LIBM: libopenlibm
LLVM: libLLVM-13.0.1 (ORCJIT, apple-m1)
Threads: 1 on 8 virtual cores
I’ve tried running the code on my computer and got exactly the same results (for the first and second run). So, things are stable and reproducible - and not random.
A third run, gives the same result as the second. So the problem should be in some initialization/compilation on the first run.
Not sure what the problem is in my system. To check - I reran it three times, and these are the exact results. Same for first two times then different in 3rd times
julia> using Turing, StableRNGs
julia> @model function example(x)
m ~ Normal(0, 1)
x ~ Normal(m, sqrt(1))
end
example (generic function with 2 methods)
julia> res1 = sample(StableRNG(123), example(1), NUTS(), 1000)
┌ Info: Found initial step size
└ ϵ = 1.6
Chains MCMC chain (1000×13×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 0.49 seconds
Compute duration = 0.49 seconds
parameters = m
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ess_per_sec
Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64
m 0.4905 0.6835 0.0216 0.0314 490.8002 1.0054 1005.7380
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
m -0.8832 0.0403 0.4865 0.9456 1.7660
julia> res2 = sample(StableRNG(123), example(1), NUTS(), 1000)
┌ Info: Found initial step size
└ ϵ = 1.6
Chains MCMC chain (1000×13×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 0.07 seconds
Compute duration = 0.07 seconds
parameters = m
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ess_per_sec
Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64
m 0.4905 0.6835 0.0216 0.0314 490.8002 1.0054 7011.4309
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
m -0.8832 0.0403 0.4865 0.9456 1.7660
julia> res3 = sample(StableRNG(123), example(1), NUTS(), 1000)
┌ Info: Found initial step size
└ ϵ = 3.2
Chains MCMC chain (1000×13×1 Array{Float64, 3}):
Iterations = 501:1:1500
Number of chains = 1
Samples per chain = 1000
Wall duration = 0.05 seconds
Compute duration = 0.05 seconds
parameters = m
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std naive_se mcse ess rhat ess_per_sec
Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64
m 0.5167 0.7029 0.0222 0.0301 533.6340 0.9999 10262.1917
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5%
Symbol Float64 Float64 Float64 Float64 Float64
m -0.8081 0.0530 0.5059 0.9911 1.8123
The problem seems to be in the initial step_size estimation. The alternative result is always given when the initial step_size is alternative as well. In any case, the following:
sample(StableRNG(123), example(1), NUTS(1000, 0.65; init_ϵ = 1.6), 1000)
sets the step size and now only one result is output.
Will investigate further.
Oh Danny Boy 
thank you so so much <3
You are welcome. BTW, found the bug:
trajectory.jl:770 in AdvancedHMC:
return find_good_stepsize(GLOBAL_RNG, h, θ; max_n_iters=max_n_iters)
uses the global RNG.
I will file an issue.
Aah - I see. Thanks a lot.