diff --git a/test/mcmc/Inference.jl b/test/mcmc/Inference.jl
index dbd14bb90..f27a9466e 100644
--- a/test/mcmc/Inference.jl
+++ b/test/mcmc/Inference.jl
@@ -19,7 +19,7 @@ using Turing
 
 @testset "Testing inference.jl with $adbackend" for adbackend in ADUtils.adbackends
     @info "Starting Inference.jl tests with $adbackend"
-    rng = StableRNG(23)
+    seed = 23
 
     @testset "threaded sampling" begin
         # Test that chains with the same seed will sample identically.
@@ -61,7 +61,7 @@ using Turing
         # Smoke test for default sample call.
         @testset "gdemo_default" begin
             chain = sample(
-                copy(rng),
+                StableRNG(seed),
                 gdemo_default,
                 HMC(0.1, 7; adtype=adbackend),
                 MCMCThreads(),
@@ -73,7 +73,7 @@ using Turing
             # run sampler: progress logging should be disabled and
             # it should return a Chains object
             sampler = Sampler(HMC(0.1, 7; adtype=adbackend), gdemo_default)
-            chains = sample(copy(rng), gdemo_default, sampler, MCMCThreads(), 10, 4)
+            chains = sample(StableRNG(seed), gdemo_default, sampler, MCMCThreads(), 10, 4)
             @test chains isa MCMCChains.Chains
         end
     end
@@ -83,29 +83,39 @@ using Turing
         alg2 = PG(20)
         alg3 = Gibbs(PG(30, :s), HMC(0.2, 4, :m; adtype=adbackend))
 
-        chn1 = sample(copy(rng), gdemo_default, alg1, 2_000; save_state=true)
+        chn1 = sample(StableRNG(seed), gdemo_default, alg1, 2_000; save_state=true)
         check_gdemo(chn1)
 
-        chn1_contd = sample(copy(rng), gdemo_default, alg1, 2_000; resume_from=chn1)
+        chn1_contd = sample(StableRNG(seed), gdemo_default, alg1, 2_000; resume_from=chn1)
         check_gdemo(chn1_contd)
 
-        chn1_contd2 = sample(copy(rng), gdemo_default, alg1, 2_000; resume_from=chn1)
+        chn1_contd2 = sample(StableRNG(seed), gdemo_default, alg1, 2_000; resume_from=chn1)
         check_gdemo(chn1_contd2)
 
         chn2 = sample(
-            copy(rng), gdemo_default, alg2, 2_000; discard_initial=100, save_state=true
+            StableRNG(seed),
+            gdemo_default,
+            alg2,
+            2_000;
+            discard_initial=100,
+            save_state=true,
         )
         check_gdemo(chn2)
 
-        chn2_contd = sample(copy(rng), gdemo_default, alg2, 2_000; resume_from=chn2)
+        chn2_contd = sample(StableRNG(seed), gdemo_default, alg2, 2_000; resume_from=chn2)
         check_gdemo(chn2_contd)
 
         chn3 = sample(
-            copy(rng), gdemo_default, alg3, 2_000; discard_initial=100, save_state=true
+            StableRNG(seed),
+            gdemo_default,
+            alg3,
+            2_000;
+            discard_initial=100,
+            save_state=true,
         )
         check_gdemo(chn3)
 
-        chn3_contd = sample(copy(rng), gdemo_default, alg3, 2_000; resume_from=chn3)
+        chn3_contd = sample(StableRNG(seed), gdemo_default, alg3, 2_000; resume_from=chn3)
         check_gdemo(chn3_contd)
     end
 
@@ -144,7 +154,7 @@ using Turing
 
         # Note that all chains contain 3 values per sample: 2 variables + log probability
         @testset "Single-threaded vanilla" begin
-            chains = sample(copy(rng), gdemo_d(), Prior(), N)
+            chains = sample(StableRNG(seed), gdemo_d(), Prior(), N)
             @test chains isa MCMCChains.Chains
             @test size(chains) == (N, 3, 1)
             @test mean(chains, :s) ≈ 3 atol = 0.1
@@ -152,7 +162,7 @@ using Turing
         end
 
         @testset "Multi-threaded" begin
-            chains = sample(copy(rng), gdemo_d(), Prior(), MCMCThreads(), N, 4)
+            chains = sample(StableRNG(seed), gdemo_d(), Prior(), MCMCThreads(), N, 4)
             @test chains isa MCMCChains.Chains
             @test size(chains) == (N, 3, 4)
             @test mean(chains, :s) ≈ 3 atol = 0.1
@@ -160,7 +170,9 @@ using Turing
         end
 
         @testset "Vector chain_type" begin
-            chains = sample(copy(rng), gdemo_d(), Prior(), N; chain_type=Vector{NamedTuple})
+            chains = sample(
+                StableRNG(seed), gdemo_d(), Prior(), N; chain_type=Vector{NamedTuple}
+            )
             @test chains isa Vector{<:NamedTuple}
             @test length(chains) == N
             @test all(length(x) == 3 for x in chains)
@@ -180,27 +192,31 @@ using Turing
             end
 
             model = issue2169_model()
-            chain = sample(copy(rng), model, Prior(), 10)
+            chain = sample(StableRNG(seed), model, Prior(), 10)
             @test all(mean(chain[:x]) .< 5)
         end
     end
 
     @testset "chain ordering" begin
         for alg in (Prior(), Emcee(10, 2.0))
-            chain_sorted = sample(copy(rng), gdemo_default, alg, 1; sort_chain=true)
+            chain_sorted = sample(StableRNG(seed), gdemo_default, alg, 1; sort_chain=true)
             @test names(MCMCChains.get_sections(chain_sorted, :parameters)) == [:m, :s]
 
-            chain_unsorted = sample(copy(rng), gdemo_default, alg, 1; sort_chain=false)
+            chain_unsorted = sample(
+                StableRNG(seed), gdemo_default, alg, 1; sort_chain=false
+            )
             @test names(MCMCChains.get_sections(chain_unsorted, :parameters)) == [:s, :m]
         end
     end
 
     @testset "chain iteration numbers" begin
         for alg in (Prior(), Emcee(10, 2.0))
-            chain = sample(copy(rng), gdemo_default, alg, 10)
+            chain = sample(StableRNG(seed), gdemo_default, alg, 10)
             @test range(chain) == 1:10
 
-            chain = sample(copy(rng), gdemo_default, alg, 10; discard_initial=5, thinning=2)
+            chain = sample(
+                StableRNG(seed), gdemo_default, alg, 10; discard_initial=5, thinning=2
+            )
             @test range(chain) == range(6; step=2, length=10)
         end
     end
@@ -216,8 +232,8 @@ using Turing
         smc = SMC()
         pg = PG(10)
 
-        res1 = sample(copy(rng), test_assume(), smc, 1_000)
-        res2 = sample(copy(rng), test_assume(), pg, 1_000)
+        res1 = sample(StableRNG(seed), test_assume(), smc, 1_000)
+        res2 = sample(StableRNG(seed), test_assume(), pg, 1_000)
 
         check_numerical(res1, [:y], [0.5]; atol=0.1)
         check_numerical(res2, [:y], [0.5]; atol=0.1)
@@ -246,9 +262,9 @@ using Turing
         pg = PG(10)
         gibbs = Gibbs(HMC(0.2, 3, :p; adtype=adbackend), PG(10, :x))
 
-        chn_s = sample(copy(rng), testbb(obs), smc, 200)
-        chn_p = sample(copy(rng), testbb(obs), pg, 200)
-        chn_g = sample(copy(rng), testbb(obs), gibbs, 200)
+        chn_s = sample(StableRNG(seed), testbb(obs), smc, 200)
+        chn_p = sample(StableRNG(seed), testbb(obs), pg, 200)
+        chn_g = sample(StableRNG(seed), testbb(obs), gibbs, 200)
 
         check_numerical(chn_s, [:p], [meanp]; atol=0.05)
         check_numerical(chn_p, [:x], [meanp]; atol=0.1)
@@ -273,7 +289,7 @@ using Turing
         end
 
         gibbs = Gibbs(PG(10, :s), HMC(0.4, 8, :m; adtype=adbackend))
-        chain = sample(copy(rng), fggibbstest(xs), gibbs, 2)
+        chain = sample(StableRNG(seed), fggibbstest(xs), gibbs, 2)
     end
 
     @testset "new grammar" begin
@@ -289,8 +305,8 @@ using Turing
             return priors
         end
 
-        chain = sample(copy(rng), gauss(x), PG(10), 10)
-        chain = sample(copy(rng), gauss(x), SMC(), 10)
+        chain = sample(StableRNG(seed), gauss(x), PG(10), 10)
+        chain = sample(StableRNG(seed), gauss(x), SMC(), 10)
 
         @model function gauss2(::Type{TV}=Vector{Float64}; x) where {TV}
             priors = TV(undef, 2)
@@ -302,14 +318,18 @@ using Turing
             return priors
         end
 
-        @test_throws ErrorException chain = sample(copy(rng), gauss2(; x=x), PG(10), 10)
-        @test_throws ErrorException chain = sample(copy(rng), gauss2(; x=x), SMC(), 10)
+        @test_throws ErrorException chain = sample(
+            StableRNG(seed), gauss2(; x=x), PG(10), 10
+        )
+        @test_throws ErrorException chain = sample(
+            StableRNG(seed), gauss2(; x=x), SMC(), 10
+        )
 
         @test_throws ErrorException chain = sample(
-            copy(rng), gauss2(DynamicPPL.TypeWrap{Vector{Float64}}(); x=x), PG(10), 10
+            StableRNG(seed), gauss2(DynamicPPL.TypeWrap{Vector{Float64}}(); x=x), PG(10), 10
         )
         @test_throws ErrorException chain = sample(
-            copy(rng), gauss2(DynamicPPL.TypeWrap{Vector{Float64}}(); x=x), SMC(), 10
+            StableRNG(seed), gauss2(DynamicPPL.TypeWrap{Vector{Float64}}(); x=x), SMC(), 10
         )
 
         @model function gauss3(x, ::Type{TV}=Vector{Float64}) where {TV}
@@ -322,13 +342,15 @@ using Turing
             return priors
         end
 
-        chain = sample(copy(rng), gauss3(x), PG(10), 10)
-        chain = sample(copy(rng), gauss3(x), SMC(), 10)
+        chain = sample(StableRNG(seed), gauss3(x), PG(10), 10)
+        chain = sample(StableRNG(seed), gauss3(x), SMC(), 10)
 
         chain = sample(
-            copy(rng), gauss3(x, DynamicPPL.TypeWrap{Vector{Real}}()), PG(10), 10
+            StableRNG(seed), gauss3(x, DynamicPPL.TypeWrap{Vector{Real}}()), PG(10), 10
+        )
+        chain = sample(
+            StableRNG(seed), gauss3(x, DynamicPPL.TypeWrap{Vector{Real}}()), SMC(), 10
         )
-        chain = sample(copy(rng), gauss3(x, DynamicPPL.TypeWrap{Vector{Real}}()), SMC(), 10)
     end
 
     # TODO(mhauru) What is this testing? Why does it not use the looped-over adbackend?
@@ -344,7 +366,7 @@ using Turing
         end
 
         sample(
-            copy(rng),
+            StableRNG(seed),
             newinterface(obs),
             HMC(0.75, 3, :p, :x; adtype=Turing.AutoForwardDiff(; chunksize=2)),
             100,
@@ -360,7 +382,9 @@ using Turing
             end
         end
 
-        chain = sample(copy(rng), noreturn([1.5 2.0]), HMC(0.1, 10; adtype=adbackend), 4000)
+        chain = sample(
+            StableRNG(seed), noreturn([1.5 2.0]), HMC(0.1, 10; adtype=adbackend), 4000
+        )
         check_numerical(chain, [:s, :m], [49 / 24, 7 / 6])
     end
 
@@ -377,9 +401,9 @@ using Turing
         smc = SMC()
         pg = PG(10)
 
-        res_is = sample(copy(rng), test(), is, 1_000)
-        res_smc = sample(copy(rng), test(), smc, 1_000)
-        res_pg = sample(copy(rng), test(), pg, 100)
+        res_is = sample(StableRNG(seed), test(), is, 1_000)
+        res_smc = sample(StableRNG(seed), test(), smc, 1_000)
+        res_pg = sample(StableRNG(seed), test(), pg, 100)
 
         @test all(isone, res_is[:x])
         @test res_is.logevidence ≈ 2 * log(0.5)
@@ -392,7 +416,7 @@ using Turing
 
     @testset "sample" begin
         alg = Gibbs(HMC(0.2, 3, :m; adtype=adbackend), PG(10, :s))
-        chn = sample(copy(rng), gdemo_default, alg, 10)
+        chn = sample(StableRNG(seed), gdemo_default, alg, 10)
     end
 
     @testset "vectorization @." begin
@@ -405,7 +429,7 @@ using Turing
 
         alg = HMC(0.01, 5; adtype=adbackend)
         x = randn(100)
-        res = sample(copy(rng), vdemo1(x), alg, 10)
+        res = sample(StableRNG(seed), vdemo1(x), alg, 10)
 
         @model function vdemo1b(x)
             s ~ InverseGamma(2, 3)
@@ -414,7 +438,7 @@ using Turing
             return s, m
         end
 
-        res = sample(copy(rng), vdemo1b(x), alg, 10)
+        res = sample(StableRNG(seed), vdemo1b(x), alg, 10)
 
         @model function vdemo2(x)
             μ ~ MvNormal(zeros(size(x, 1)), I)
@@ -423,7 +447,7 @@ using Turing
 
         D = 2
         alg = HMC(0.01, 5; adtype=adbackend)
-        res = sample(copy(rng), vdemo2(randn(D, 100)), alg, 10)
+        res = sample(StableRNG(seed), vdemo2(randn(D, 100)), alg, 10)
 
         # Vector assumptions
         N = 10
@@ -438,7 +462,7 @@ using Turing
 
         # TODO(mhauru) What is the point of the below @elapsed stuff? It prints out some
         # timings. Do we actually ever look at them?
-        t_loop = @elapsed res = sample(copy(rng), vdemo3(), alg, 1000)
+        t_loop = @elapsed res = sample(StableRNG(seed), vdemo3(), alg, 1000)
 
         # Test for vectorize UnivariateDistribution
         @model function vdemo4()
@@ -446,11 +470,11 @@ using Turing
             @. x ~ Normal(0, 2)
         end
 
-        t_vec = @elapsed res = sample(copy(rng), vdemo4(), alg, 1000)
+        t_vec = @elapsed res = sample(StableRNG(seed), vdemo4(), alg, 1000)
 
         @model vdemo5() = x ~ MvNormal(zeros(N), 4 * I)
 
-        t_mv = @elapsed res = sample(copy(rng), vdemo5(), alg, 1000)
+        t_mv = @elapsed res = sample(StableRNG(seed), vdemo5(), alg, 1000)
 
         println("Time for")
         println("  Loop : ", t_loop)
@@ -463,7 +487,7 @@ using Turing
             @. x ~ InverseGamma(2, 3)
         end
 
-        sample(copy(rng), vdemo6(), alg, 10)
+        sample(StableRNG(seed), vdemo6(), alg, 10)
 
         N = 3
         @model function vdemo7()
@@ -471,7 +495,7 @@ using Turing
             @. x ~ [InverseGamma(2, 3) for i in 1:N]
         end
 
-        sample(copy(rng), vdemo7(), alg, 10)
+        sample(StableRNG(seed), vdemo7(), alg, 10)
     end
 
     @testset "vectorization .~" begin
@@ -484,7 +508,7 @@ using Turing
 
         alg = HMC(0.01, 5; adtype=adbackend)
         x = randn(100)
-        res = sample(copy(rng), vdemo1(x), alg, 10)
+        res = sample(StableRNG(seed), vdemo1(x), alg, 10)
 
         @model function vdemo2(x)
             μ ~ MvNormal(zeros(size(x, 1)), I)
@@ -493,7 +517,7 @@ using Turing
 
         D = 2
         alg = HMC(0.01, 5; adtype=adbackend)
-        res = sample(copy(rng), vdemo2(randn(D, 100)), alg, 10)
+        res = sample(StableRNG(seed), vdemo2(randn(D, 100)), alg, 10)
 
         # Vector assumptions
         N = 10
@@ -507,7 +531,7 @@ using Turing
         end
 
         # TODO(mhauru) Same question as above about @elapsed.
-        t_loop = @elapsed res = sample(copy(rng), vdemo3(), alg, 1_000)
+        t_loop = @elapsed res = sample(StableRNG(seed), vdemo3(), alg, 1_000)
 
         # Test for vectorize UnivariateDistribution
         @model function vdemo4()
@@ -515,11 +539,11 @@ using Turing
             return x .~ Normal(0, 2)
         end
 
-        t_vec = @elapsed res = sample(copy(rng), vdemo4(), alg, 1_000)
+        t_vec = @elapsed res = sample(StableRNG(seed), vdemo4(), alg, 1_000)
 
         @model vdemo5() = x ~ MvNormal(zeros(N), 4 * I)
 
-        t_mv = @elapsed res = sample(copy(rng), vdemo5(), alg, 1_000)
+        t_mv = @elapsed res = sample(StableRNG(seed), vdemo5(), alg, 1_000)
 
         println("Time for")
         println("  Loop : ", t_loop)
@@ -532,14 +556,14 @@ using Turing
             return x .~ InverseGamma(2, 3)
         end
 
-        sample(copy(rng), vdemo6(), alg, 10)
+        sample(StableRNG(seed), vdemo6(), alg, 10)
 
         @model function vdemo7()
             x = Array{Real}(undef, N, N)
             return x .~ [InverseGamma(2, 3) for i in 1:N]
         end
 
-        sample(copy(rng), vdemo7(), alg, 10)
+        sample(StableRNG(seed), vdemo7(), alg, 10)
     end
 
     @testset "Type parameters" begin
@@ -555,33 +579,35 @@ using Turing
 
         # TODO(mhauru) What are we testing below? Just that using a type parameter doesn't
         # crash?
-        sample(copy(rng), vdemo1(), alg, 10)
-        sample(copy(rng), vdemo1(DynamicPPL.TypeWrap{Float64}()), alg, 10)
+        sample(StableRNG(seed), vdemo1(), alg, 10)
+        sample(StableRNG(seed), vdemo1(DynamicPPL.TypeWrap{Float64}()), alg, 10)
 
         vdemo1kw(; T) = vdemo1(T)
-        sample(copy(rng), vdemo1kw(; T=DynamicPPL.TypeWrap{Float64}()), alg, 10)
+        sample(StableRNG(seed), vdemo1kw(; T=DynamicPPL.TypeWrap{Float64}()), alg, 10)
 
         @model function vdemo2(::Type{T}=Float64) where {T<:Real}
             x = Vector{T}(undef, N)
             @. x ~ Normal(0, 2)
         end
 
-        sample(copy(rng), vdemo2(), alg, 10)
-        sample(copy(rng), vdemo2(DynamicPPL.TypeWrap{Float64}()), alg, 10)
+        sample(StableRNG(seed), vdemo2(), alg, 10)
+        sample(StableRNG(seed), vdemo2(DynamicPPL.TypeWrap{Float64}()), alg, 10)
 
         vdemo2kw(; T) = vdemo2(T)
-        sample(copy(rng), vdemo2kw(; T=DynamicPPL.TypeWrap{Float64}()), alg, 10)
+        sample(StableRNG(seed), vdemo2kw(; T=DynamicPPL.TypeWrap{Float64}()), alg, 10)
 
         @model function vdemo3(::Type{TV}=Vector{Float64}) where {TV<:AbstractVector}
             x = TV(undef, N)
             @. x ~ InverseGamma(2, 3)
         end
 
-        sample(copy(rng), vdemo3(), alg, 10)
-        sample(copy(rng), vdemo3(DynamicPPL.TypeWrap{Vector{Float64}}()), alg, 10)
+        sample(StableRNG(seed), vdemo3(), alg, 10)
+        sample(StableRNG(seed), vdemo3(DynamicPPL.TypeWrap{Vector{Float64}}()), alg, 10)
 
         vdemo3kw(; T) = vdemo3(T)
-        sample(copy(rng), vdemo3kw(; T=DynamicPPL.TypeWrap{Vector{Float64}}()), alg, 10)
+        sample(
+            StableRNG(seed), vdemo3kw(; T=DynamicPPL.TypeWrap{Vector{Float64}}()), alg, 10
+        )
     end
 
     @testset "names_values" begin
@@ -597,17 +623,21 @@ using Turing
         end
 
         @test_throws ErrorException sample(
-            copy(rng), demo_repeated_varname(), NUTS(), 10; check_model=true
+            StableRNG(seed), demo_repeated_varname(), NUTS(), 10; check_model=true
         )
         # Make sure that disabling the check also works.
-        @test (sample(copy(rng), demo_repeated_varname(), Prior(), 10; check_model=false);
-        true)
+        @test (
+            sample(
+                StableRNG(seed), demo_repeated_varname(), Prior(), 10; check_model=false
+            );
+            true
+        )
 
         @model function demo_incorrect_missing(y)
             return y[1:1] ~ MvNormal(zeros(1), I)
         end
         @test_throws ErrorException sample(
-            copy(rng), demo_incorrect_missing([missing]), NUTS(), 10; check_model=true
+            StableRNG(seed), demo_incorrect_missing([missing]), NUTS(), 10; check_model=true
         )
     end
 end
diff --git a/test/mcmc/ess.jl b/test/mcmc/ess.jl
index edd0b9959..23f4a11ae 100644
--- a/test/mcmc/ess.jl
+++ b/test/mcmc/ess.jl
@@ -46,27 +46,27 @@ using Turing
 
     @testset "ESS inference" begin
         @info "Starting ESS inference tests"
-        rng = StableRNG(23)
+        seed = 23
 
         @testset "demo_default" begin
-            chain = sample(copy(rng), demo_default, ESS(), 5_000)
+            chain = sample(StableRNG(seed), demo_default, ESS(), 5_000)
             check_numerical(chain, [:m], [0.8]; atol=0.1)
         end
 
         @testset "demodot_default" begin
-            chain = sample(copy(rng), demodot_default, ESS(), 5_000)
+            chain = sample(StableRNG(seed), demodot_default, ESS(), 5_000)
             check_numerical(chain, ["m[1]", "m[2]"], [0.0, 0.8]; atol=0.1)
         end
 
         @testset "gdemo with CSMC + ESS" begin
             alg = Gibbs(CSMC(15, :s), ESS(:m))
-            chain = sample(copy(rng), gdemo(1.5, 2.0), alg, 2000)
+            chain = sample(StableRNG(seed), gdemo(1.5, 2.0), alg, 2000)
             check_numerical(chain, [:s, :m], [49 / 24, 7 / 6]; atol=0.1)
         end
 
         @testset "MoGtest_default with CSMC + ESS" begin
             alg = Gibbs(CSMC(15, :z1, :z2, :z3, :z4), ESS(:mu1), ESS(:mu2))
-            chain = sample(copy(rng), MoGtest_default, alg, 2000)
+            chain = sample(StableRNG(seed), MoGtest_default, alg, 2000)
             check_MoGtest_default(chain; atol=0.1)
         end
 
diff --git a/test/mcmc/hmc.jl b/test/mcmc/hmc.jl
index c25104063..7d3f3a380 100644
--- a/test/mcmc/hmc.jl
+++ b/test/mcmc/hmc.jl
@@ -20,7 +20,7 @@ using Turing
 
 @testset "Testing hmc.jl with $adbackend" for adbackend in ADUtils.adbackends
     @info "Starting HMC tests with $adbackend"
-    rng = StableRNG(123)
+    seed = 123
 
     @testset "constrained bounded" begin
         obs = [0, 1, 0, 1, 1, 1, 1, 1, 1, 1]
@@ -34,7 +34,7 @@ using Turing
         end
 
         chain = sample(
-            copy(rng),
+            StableRNG(seed),
             constrained_test(obs),
             HMC(1.5, 3; adtype=adbackend),# using a large step size (1.5)
             1_000,
@@ -56,7 +56,10 @@ using Turing
         end
 
         chain = sample(
-            copy(rng), constrained_simplex_test(obs12), HMC(0.75, 2; adtype=adbackend), 1000
+            StableRNG(seed),
+            constrained_simplex_test(obs12),
+            HMC(0.75, 2; adtype=adbackend),
+            1000,
         )
 
         check_numerical(chain, ["ps[1]", "ps[2]"], [5 / 16, 11 / 16]; atol=0.015)
@@ -64,7 +67,7 @@ using Turing
 
     @testset "hmc reverse diff" begin
         alg = HMC(0.1, 10; adtype=adbackend)
-        res = sample(copy(rng), gdemo_default, alg, 4_000)
+        res = sample(StableRNG(seed), gdemo_default, alg, 4_000)
         check_gdemo(res; rtol=0.1)
     end
 
@@ -129,12 +132,12 @@ using Turing
         end
 
         # Sampling
-        chain = sample(copy(rng), bnn(ts), HMC(0.1, 5; adtype=adbackend), 10)
+        chain = sample(StableRNG(seed), bnn(ts), HMC(0.1, 5; adtype=adbackend), 10)
     end
 
     @testset "hmcda inference" begin
         alg1 = HMCDA(500, 0.8, 0.015; adtype=adbackend)
-        res1 = sample(copy(rng), gdemo_default, alg1, 3_000)
+        res1 = sample(StableRNG(seed), gdemo_default, alg1, 3_000)
         check_gdemo(res1)
     end
 
@@ -167,7 +170,7 @@ using Turing
 
     @testset "nuts inference" begin
         alg = NUTS(1000, 0.8; adtype=adbackend)
-        res = sample(copy(rng), gdemo_default, alg, 500)
+        res = sample(StableRNG(seed), gdemo_default, alg, 500)
         check_gdemo(res)
     end
 
@@ -188,8 +191,8 @@ using Turing
     @testset "check discard" begin
         alg = NUTS(100, 0.8; adtype=adbackend)
 
-        c1 = sample(copy(rng), gdemo_default, alg, 500; discard_adapt=true)
-        c2 = sample(copy(rng), gdemo_default, alg, 500; discard_adapt=false)
+        c1 = sample(StableRNG(seed), gdemo_default, alg, 500; discard_adapt=true)
+        c2 = sample(StableRNG(seed), gdemo_default, alg, 500; discard_adapt=false)
 
         @test size(c1, 1) == 500
         @test size(c2, 1) == 500
@@ -199,9 +202,9 @@ using Turing
         alg1 = Gibbs(PG(10, :m), NUTS(100, 0.65, :s; adtype=adbackend))
         alg2 = Gibbs(PG(10, :m), HMC(0.1, 3, :s; adtype=adbackend))
         alg3 = Gibbs(PG(10, :m), HMCDA(100, 0.65, 0.3, :s; adtype=adbackend))
-        @test sample(copy(rng), gdemo_default, alg1, 10) isa Chains
-        @test sample(copy(rng), gdemo_default, alg2, 10) isa Chains
-        @test sample(copy(rng), gdemo_default, alg3, 10) isa Chains
+        @test sample(StableRNG(seed), gdemo_default, alg1, 10) isa Chains
+        @test sample(StableRNG(seed), gdemo_default, alg2, 10) isa Chains
+        @test sample(StableRNG(seed), gdemo_default, alg3, 10) isa Chains
     end
 
     @testset "Regression tests" begin
@@ -210,28 +213,28 @@ using Turing
             m = Matrix{T}(undef, 2, 3)
             return m .~ MvNormal(zeros(2), I)
         end
-        @test sample(copy(rng), mwe1(), HMC(0.2, 4; adtype=adbackend), 100) isa Chains
+        @test sample(StableRNG(seed), mwe1(), HMC(0.2, 4; adtype=adbackend), 100) isa Chains
 
         @model function mwe2(::Type{T}=Matrix{Float64}) where {T}
             m = T(undef, 2, 3)
             return m .~ MvNormal(zeros(2), I)
         end
-        @test sample(copy(rng), mwe2(), HMC(0.2, 4; adtype=adbackend), 100) isa Chains
+        @test sample(StableRNG(seed), mwe2(), HMC(0.2, 4; adtype=adbackend), 100) isa Chains
 
         # https://github.com/TuringLang/Turing.jl/issues/1308
         @model function mwe3(::Type{T}=Array{Float64}) where {T}
             m = T(undef, 2, 3)
             return m .~ MvNormal(zeros(2), I)
         end
-        @test sample(copy(rng), mwe3(), HMC(0.2, 4; adtype=adbackend), 100) isa Chains
+        @test sample(StableRNG(seed), mwe3(), HMC(0.2, 4; adtype=adbackend), 100) isa Chains
     end
 
     # issue #1923
     @testset "reproducibility" begin
         alg = NUTS(1000, 0.8; adtype=adbackend)
-        res1 = sample(copy(rng), gdemo_default, alg, 10)
-        res2 = sample(copy(rng), gdemo_default, alg, 10)
-        res3 = sample(copy(rng), gdemo_default, alg, 10)
+        res1 = sample(StableRNG(seed), gdemo_default, alg, 10)
+        res2 = sample(StableRNG(seed), gdemo_default, alg, 10)
+        res3 = sample(StableRNG(seed), gdemo_default, alg, 10)
         @test Array(res1) == Array(res2) == Array(res3)
     end
 
@@ -332,7 +335,7 @@ using Turing
             gdemo_default, ADTypeCheckContext(adbackend, gdemo_default.context)
         )
         # These will error if the adbackend being used is not the one set.
-        sample(copy(rng), m, alg, 10)
+        sample(StableRNG(seed), m, alg, 10)
     end
 end
 
diff --git a/test/mcmc/mh.jl b/test/mcmc/mh.jl
index 191d9646e..d71a5fbc6 100644
--- a/test/mcmc/mh.jl
+++ b/test/mcmc/mh.jl
@@ -19,7 +19,7 @@ GKernel(var) = (x) -> Normal(x, sqrt.(var))
 
 @testset "mh.jl" begin
     @info "Starting MH tests"
-    rng = StableRNG(23)
+    seed = 23
 
     @testset "mh constructor" begin
         N = 10
@@ -55,7 +55,7 @@ GKernel(var) = (x) -> Normal(x, sqrt.(var))
         @testset "gdemo_default" begin
             alg = MH()
             chain = sample(
-                copy(rng), gdemo_default, alg, 10_000; discard_initial, initial_params
+                StableRNG(seed), gdemo_default, alg, 10_000; discard_initial, initial_params
             )
             check_gdemo(chain; atol=0.1)
         end
@@ -63,7 +63,7 @@ GKernel(var) = (x) -> Normal(x, sqrt.(var))
         @testset "gdemo_default with custom proposals" begin
             alg = MH((:s, InverseGamma(2, 3)), (:m, GKernel(1.0)))
             chain = sample(
-                copy(rng), gdemo_default, alg, 10_000; discard_initial, initial_params
+                StableRNG(seed), gdemo_default, alg, 10_000; discard_initial, initial_params
             )
             check_gdemo(chain; atol=0.1)
         end
@@ -71,7 +71,7 @@ GKernel(var) = (x) -> Normal(x, sqrt.(var))
         @testset "gdemo_default with MH-within-Gibbs" begin
             alg = Gibbs(MH(:m), MH(:s))
             chain = sample(
-                copy(rng), gdemo_default, alg, 10_000; discard_initial, initial_params
+                StableRNG(seed), gdemo_default, alg, 10_000; discard_initial, initial_params
             )
             check_gdemo(chain; atol=0.1)
         end
@@ -81,7 +81,7 @@ GKernel(var) = (x) -> Normal(x, sqrt.(var))
                 CSMC(15, :z1, :z2, :z3, :z4), MH((:mu1, GKernel(1)), (:mu2, GKernel(1)))
             )
             chain = sample(
-                copy(rng),
+                StableRNG(seed),
                 MoGtest_default,
                 gibbs,
                 500;
@@ -145,8 +145,8 @@ GKernel(var) = (x) -> Normal(x, sqrt.(var))
         @test spl1.proposals.proposal.Σ.mat == spl2.proposals.proposal.Σ.mat
 
         # Test inference.
-        chain1 = sample(copy(rng), gdemo_default, spl1, 2_000)
-        chain2 = sample(copy(rng), gdemo_default, spl2, 2_000)
+        chain1 = sample(StableRNG(seed), gdemo_default, spl1, 2_000)
+        chain2 = sample(StableRNG(seed), gdemo_default, spl2, 2_000)
 
         check_gdemo(chain1)
         check_gdemo(chain2)
@@ -178,8 +178,8 @@ GKernel(var) = (x) -> Normal(x, sqrt.(var))
         alg_small = Gibbs(MH((:μ, vc_μ)), MH((:σ, vc_σ)))
         alg_big = MH()
 
-        chn_small = sample(copy(rng), mod, alg_small, 1_000)
-        chn_big = sample(copy(rng), mod, alg_big, 1_000)
+        chn_small = sample(StableRNG(seed), mod, alg_small, 1_000)
+        chn_big = sample(StableRNG(seed), mod, alg_big, 1_000)
 
         # Test that the small variance version is actually smaller.
         variance_small = var(diff(Array(chn_small["μ[1]"]); dims=1))
@@ -195,12 +195,12 @@ GKernel(var) = (x) -> Normal(x, sqrt.(var))
             end
         end
 
-        chain = sample(copy(rng), test(1), MH(), 5_000)
+        chain = sample(StableRNG(seed), test(1), MH(), 5_000)
         for i in 1:5
             @test mean(chain, "T[1][$i]") ≈ 0.2 atol = 0.01
         end
 
-        chain = sample(copy(rng), test(10), MH(), 5_000)
+        chain = sample(StableRNG(seed), test(10), MH(), 5_000)
         for j in 1:10, i in 1:5
             @test mean(chain, "T[$j][$i]") ≈ 0.2 atol = 0.01
         end
@@ -209,7 +209,7 @@ GKernel(var) = (x) -> Normal(x, sqrt.(var))
     @testset "LKJCholesky" begin
         for uplo in ['L', 'U']
             @model f() = x ~ LKJCholesky(2, 1, uplo)
-            chain = sample(copy(rng), f(), MH(), 5_000)
+            chain = sample(StableRNG(seed), f(), MH(), 5_000)
             indices = [(1, 1), (2, 1), (2, 2)]
             values = [1, 0, 0.785]
             for ((i, j), v) in zip(indices, values)
@@ -270,7 +270,12 @@ GKernel(var) = (x) -> Normal(x, sqrt.(var))
         burnin = 10_000
         n = 10_000
         chain = sample(
-            copy(rng), gdemo_default_prior, alg, n; discard_initial=burnin, thinning=10
+            StableRNG(seed),
+            gdemo_default_prior,
+            alg,
+            n;
+            discard_initial=burnin,
+            thinning=10,
         )
         check_numerical(chain, [:s, :m], [mean(InverseGamma(2, 3)), 0]; atol=0.3)
     end
@@ -278,7 +283,7 @@ GKernel(var) = (x) -> Normal(x, sqrt.(var))
     @testset "`filldist` proposal (issue #2180)" begin
         @model demo_filldist_issue2180() = x ~ MvNormal(zeros(3), I)
         chain = sample(
-            copy(rng),
+            StableRNG(seed),
             demo_filldist_issue2180(),
             MH(AdvancedMH.RandomWalkProposal(filldist(Normal(), 3))),
             10_000,