# StableRNG in Turing is not producing reproducible output
**URL:** https://discourse.julialang.org/t/stablerng-in-turing-is-not-producing-reproducible-output/92032
**Category:** Statistics
**Tags:** turing, reproducibility, stablerng
**Created:** [December 23, 2022, 4:41am UTC](https://discourse.julialang.org/t/stablerng-in-turing-is-not-producing-reproducible-output/92032 "2022-12-23T04:41:06Z")
**Posts on this page:** 13
**Page:** 1
### Author:  [@Sourish](https://discourse.julialang.org/u/Sourish)
#### Post date: [December 23, 2022, 4:41am UTC](https://discourse.julialang.org/t/stablerng-in-turing-is-not-producing-reproducible-output/92032/1 "2022-12-23T04:41:07Z")
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.
```julia
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)
```
---
### Author:  [@Sukera](https://discourse.julialang.org/u/Sukera)
#### Post date: [December 23, 2022, 11:04am UTC](https://discourse.julialang.org/t/stablerng-in-turing-is-not-producing-reproducible-output/92032/2 "2022-12-23T11:04:45Z")
> [@Sourish](#):
>
> However, the results are not reproducing.
In what sense? What are you changing between runs?
---
### Author:  [@ericphanson](https://discourse.julialang.org/u/ericphanson)
#### Post date: [December 23, 2022, 12:19pm UTC](https://discourse.julialang.org/t/stablerng-in-turing-is-not-producing-reproducible-output/92032/3 "2022-12-23T12:19:09Z")
I’m not really sure what to look for, but I tried running that code twice and the printed output looked the same:
```julia
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
```
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [December 23, 2022, 12:30pm UTC](https://discourse.julialang.org/t/stablerng-in-turing-is-not-producing-reproducible-output/92032/4 "2022-12-23T12:30:50Z")
I would move to single-threaded as a first step to isolate non-reproducibility.
---
### Author:  [@Sourish](https://discourse.julialang.org/u/Sourish)
#### Post date: [December 23, 2022, 2:25pm UTC](https://discourse.julialang.org/t/stablerng-in-turing-is-not-producing-reproducible-output/92032/5 "2022-12-23T14:25:18Z")
Thank you for showing interest in my queries. This is exactly what I did. Please see below.
```julia
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.
---
### Author:  [@Sourish](https://discourse.julialang.org/u/Sourish)
#### Post date: [December 23, 2022, 2:28pm UTC](https://discourse.julialang.org/t/stablerng-in-turing-is-not-producing-reproducible-output/92032/6 "2022-12-23T14:28:16Z")
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
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
```
---
### Author:  [@Sourish](https://discourse.julialang.org/u/Sourish)
#### Post date: [December 23, 2022, 2:29pm UTC](https://discourse.julialang.org/t/stablerng-in-turing-is-not-producing-reproducible-output/92032/7 "2022-12-23T14:29:46Z")
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
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
```
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [December 23, 2022, 2:51pm UTC](https://discourse.julialang.org/t/stablerng-in-turing-is-not-producing-reproducible-output/92032/8 "2022-12-23T14:51:29Z")
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.
---
### Author:  [@Sourish](https://discourse.julialang.org/u/Sourish)
#### Post date: [December 23, 2022, 3:03pm UTC](https://discourse.julialang.org/t/stablerng-in-turing-is-not-producing-reproducible-output/92032/9 "2022-12-23T15:03:49Z")
> [@ericphanson](#):
>
> `st StableRNGs Turing`
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
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
```
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [December 23, 2022, 3:18pm UTC](https://discourse.julialang.org/t/stablerng-in-turing-is-not-producing-reproducible-output/92032/10 "2022-12-23T15:18:25Z")
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:
```julia
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.
---
### Author:  [@Sourish](https://discourse.julialang.org/u/Sourish)
#### Post date: [December 23, 2022, 3:25pm UTC](https://discourse.julialang.org/t/stablerng-in-turing-is-not-producing-reproducible-output/92032/11 "2022-12-23T15:25:09Z")
> [@Dan](#):
>
> ```julia
> sample(StableRNG(123), example(1), NUTS(1000, 0.65; init_ϵ = 1.6), 1000)
>
> ```
Oh Danny Boy 😃 🙂 thank you so so much \<3
---
### Author:  [@Dan](https://discourse.julialang.org/u/Dan)
#### Post date: [December 23, 2022, 3:32pm UTC](https://discourse.julialang.org/t/stablerng-in-turing-is-not-producing-reproducible-output/92032/12 "2022-12-23T15:32:32Z")
You are welcome. BTW, found the bug:
trajectory.jl:770 in AdvancedHMC:
```julia
return find_good_stepsize(GLOBAL_RNG, h, θ; max_n_iters=max_n_iters)
```
uses the global RNG.
I will file an issue.
---
### Author:  [@Sourish](https://discourse.julialang.org/u/Sourish)
#### Post date: [December 23, 2022, 3:45pm UTC](https://discourse.julialang.org/t/stablerng-in-turing-is-not-producing-reproducible-output/92032/13 "2022-12-23T15:45:36Z")
Aah - I see. Thanks a lot.