diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index fba1f2e2..6ecac852 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -49,3 +49,14 @@ jobs:
         with:
           github-token: ${{ secrets.GITHUB_TOKEN }}
           path-to-lcov: lcov.info
+          parallel: true
+
+  finish:
+    needs: test
+    if: ${{ always() }}
+    runs-on: ubuntu-latest
+    steps:
+    - name: Coveralls Finished
+      uses: coverallsapp/github-action@v2
+      with:
+        parallel-finished: true
diff --git a/Project.toml b/Project.toml
index 68f14dfa..a046140d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -29,5 +29,5 @@ StatsBase = "0.33 - 0.34"
 StatsModels = "0.6.28 - 0.7"
 TableOperations = "1.2"
 Tables = "1.6"
-Turing = "0.20 - 0.39"
+Turing = "0.20 - 0.40"
 julia = "1"
diff --git a/test/Project.toml b/test/Project.toml
index 80eed6c8..5612711e 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -6,6 +6,7 @@ ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
diff --git a/test/runtests.jl b/test/runtests.jl
index f3d0cbfa..ddc71043 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -6,7 +6,6 @@ using CategoricalArrays: CategoricalValue
 using CategoricalArrays: categorical, levels
 using Statistics: mean, std
 using TimerOutputs: TimerOutputs, @timeit
-using Random: seed!
 
 const T = TuringGLM
 
diff --git a/test/turing_model.jl b/test/turing_model.jl
index fa833fea..1f4413d3 100644
--- a/test/turing_model.jl
+++ b/test/turing_model.jl
@@ -1,3 +1,5 @@
+using StableRNGs: StableRNG
+
 @timed_testset "turing_model" begin
     DATA_DIR = joinpath("..", "data")
     kidiq = CSV.read(joinpath(DATA_DIR, "kidiq.csv"), DataFrame)
@@ -8,7 +10,7 @@
         f = @formula(kid_score ~ mom_iq * mom_hs)
         @testset "standardize=false" begin
             m = turing_model(f, kidiq)
-            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
+            chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2)
             @test summarystats(chn)[:α, :mean] ≈ 31.80 atol = 2.0
             @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.507 atol = 0.2
             @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.22 atol = 0.2
@@ -16,7 +18,7 @@
 
         @testset "standardize=true" begin
             m = turing_model(f, kidiq; standardize=true)
-            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
+            chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2)
             @test summarystats(chn)[:α, :mean] ≈ 0.000 atol = 0.2
             @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.648 atol = 0.2
             @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.849 atol = 0.2
@@ -25,14 +27,14 @@
         @testset "custom_priors" begin
             priors = CustomPrior(Normal(), Normal(28, 5), nothing)
             m = turing_model(f, kidiq; priors)
-            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
+            chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2)
             @test summarystats(chn)[:α, :mean] ≈ 28.758 atol = 2.0
             @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.539 atol = 0.2
             @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.3863 atol = 0.2
         end
         @testset "explicit calling Normal" begin
             m = turing_model(f, kidiq; model=Normal)
-            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
+            chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2)
             @test summarystats(chn)[:α, :mean] ≈ 31.80 atol = 2.0
             @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.507 atol = 0.2
             @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.22 atol = 0.2
@@ -42,7 +44,7 @@
         f = @formula(kid_score ~ mom_iq * mom_hs)
         @testset "standardize=false" begin
             m = turing_model(f, kidiq; model=TDist)
-            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
+            chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2)
             @test summarystats(chn)[:α, :mean] ≈ 33.31 atol = 2.0
             @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.519 atol = 0.2
             @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.340 atol = 0.2
@@ -52,7 +54,7 @@
         @testset "custom_priors" begin
             priors = CustomPrior(Normal(), Normal(28, 5), Exponential(2))
             m = turing_model(f, kidiq; model=TDist, priors)
-            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
+            chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2)
             @test summarystats(chn)[:α, :mean] ≈ 28.565 atol = 2.0
             @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.551 atol = 0.2
             @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ 0.255 atol = 0.2
@@ -63,7 +65,7 @@
         f = @formula(switch ~ arsenic + dist + assoc + educ)
         @testset "standardize=false" begin
             m = turing_model(f, wells; model=Bernoulli)
-            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
+            chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2)
             @test summarystats(chn)[:α, :mean] ≈ -0.153 atol = 0.2
             @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.467 atol = 0.2
             @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.009 atol = 0.2
@@ -72,7 +74,7 @@
         @testset "custom_priors" begin
             priors = CustomPrior(Normal(), Normal(), nothing)
             m = turing_model(f, wells; model=Bernoulli, priors)
-            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
+            chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2)
             @test summarystats(chn)[:α, :mean] ≈ -0.155 atol = 0.2
             @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.468 atol = 0.2
             @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.009 atol = 0.2
@@ -82,7 +84,8 @@
         f = @formula(y ~ roach1 + treatment + senior + exposure2)
         @testset "standardize=false" begin
             m = turing_model(f, roaches; model=Poisson)
-            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
+            # seed of 123 gives bad results
+            chn = sample(StableRNG(124), m, NUTS(), MCMCThreads(), 2_000, 2)
             @test summarystats(chn)[:α, :mean] ≈ 2.969 atol = 0.5
             @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.006 atol = 0.2
             @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.5145 atol = 0.2
@@ -91,7 +94,7 @@
         @testset "custom_priors" begin
             priors = CustomPrior(Normal(2, 5), Normal(), nothing)
             m = turing_model(f, roaches; model=Poisson, priors)
-            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
+            chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2)
             @test summarystats(chn)[:α, :mean] ≈ 2.963 atol = 0.5
             @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.006 atol = 0.2
             @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.5145 atol = 0.2
@@ -101,7 +104,7 @@
         f = @formula(y ~ roach1 + treatment + senior + exposure2)
         @testset "standardize=false" begin
             m = turing_model(f, roaches; model=NegativeBinomial)
-            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
+            chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2)
             @test summarystats(chn)[:α, :mean] ≈ 2.448 atol = 0.5
             @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.013 atol = 0.2
             @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.734 atol = 0.2
@@ -111,7 +114,7 @@
         @testset "custom_priors" begin
             priors = CustomPrior(Normal(), Normal(2, 5), Exponential(0.5))
             m = turing_model(f, roaches; model=NegativeBinomial, priors)
-            chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
+            chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2)
             @test summarystats(chn)[:α, :mean] ≈ 2.401 atol = 0.5
             @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 0.013 atol = 0.2
             @test quantile(chn)[Symbol("β[2]"), Symbol("50.0%")] ≈ -0.723 atol = 0.2
@@ -121,7 +124,7 @@
     @timed_testset "Hierarchical Model" begin
         f = @formula(y ~ (1 | cheese) + background)
         m = turing_model(f, cheese)
-        chn = sample(seed!(123), m, NUTS(), MCMCThreads(), 2_000, 2)
+        chn = sample(StableRNG(123), m, NUTS(), MCMCThreads(), 2_000, 2)
         @test summarystats(chn)[:α, :mean] ≈ 68.07 atol = 2.0
         @test summarystats(chn)[Symbol("β[1]"), :mean] ≈ 6.60 atol = 0.2
         @test summarystats(chn)[Symbol("zⱼ[1]"), :mean] ≈ 0.348 atol = 0.2