diff --git a/test/Project.toml b/test/Project.toml
index c1af7b70..e931345c 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -11,6 +11,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PairPlots = "43a3c2be-4208-490b-832a-a21dcd55d7da"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
diff --git a/test/data/eigenanalysis-1.png b/test/data/eigenanalysis-1.png
index cc269614..0cc989a6 100644
Binary files a/test/data/eigenanalysis-1.png and b/test/data/eigenanalysis-1.png differ
diff --git a/test/data/eigenanalysis-2.png b/test/data/eigenanalysis-2.png
index a4755cb9..533b5ed2 100644
Binary files a/test/data/eigenanalysis-2.png and b/test/data/eigenanalysis-2.png differ
diff --git a/test/data/projectionpursuit-1.png b/test/data/projectionpursuit-1.png
index 927de853..d2ecd696 100644
Binary files a/test/data/projectionpursuit-1.png and b/test/data/projectionpursuit-1.png differ
diff --git a/test/data/projectionpursuit-2.png b/test/data/projectionpursuit-2.png
index 413b15af..fcb44516 100644
Binary files a/test/data/projectionpursuit-2.png and b/test/data/projectionpursuit-2.png differ
diff --git a/test/data/projectionpursuit-3.png b/test/data/projectionpursuit-3.png
index dc8d34e4..b21ecf0a 100644
Binary files a/test/data/projectionpursuit-3.png and b/test/data/projectionpursuit-3.png differ
diff --git a/test/runtests.jl b/test/runtests.jl
index 42f1c03d..07f06e60 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,9 +10,9 @@ using StatsBase
 using Statistics
 using DelimitedFiles
 using ReferenceTests
+using StableRNGs
 using PairPlots
 using ImageIO
-using Random
 using Test
 
 import CairoMakie as Mke
@@ -30,9 +30,8 @@ islinux = Sys.islinux()
 visualtests = !isCI || (isCI && islinux)
 datadir = joinpath(@__DIR__, "data")
 
-# using MersenneTwister for backward
-# compatibility with old Julia versions
-rng = MersenneTwister(42)
+# using StableRNG for compatibility between Julia versions
+rng = StableRNG(42)
 
 # for functor tests in Functional testset
 struct Polynomial{T<:Real}
diff --git a/test/shows.jl b/test/shows.jl
index cdd810e1..06a05f39 100644
--- a/test/shows.jl
+++ b/test/shows.jl
@@ -143,7 +143,7 @@
     ├─ weights: nothing
     ├─ replace: false
     ├─ ordered: true
-    └─ rng: TaskLocalRNG()"""
+    └─ rng: Random.TaskLocalRNG()"""
   end
 
   @testset "Filter" begin
diff --git a/test/transforms/eigenanalysis.jl b/test/transforms/eigenanalysis.jl
index f180ccb7..d2367abd 100644
--- a/test/transforms/eigenanalysis.jl
+++ b/test/transforms/eigenanalysis.jl
@@ -2,8 +2,8 @@
   @test TT.parameters(EigenAnalysis(:V)) == (proj=:V, maxdim=nothing, pratio=1.0)
 
   # PCA test
-  x = rand(Normal(0, 10), 1500)
-  y = x + rand(Normal(0, 2), 1500)
+  x = rand(rng, Normal(0, 10), 1500)
+  y = x + rand(rng, Normal(0, 2), 1500)
   t = Table(; x, y)
   T = EigenAnalysis(:V)
   n, c = apply(T, t)
@@ -16,8 +16,8 @@
   @test Tables.matrix(t) ≈ Tables.matrix(tₒ)
 
   # DRS test
-  x = rand(Normal(0, 10), 1500)
-  y = x + rand(Normal(0, 2), 1500)
+  x = rand(rng, Normal(0, 10), 1500)
+  y = x + rand(rng, Normal(0, 2), 1500)
   t = Table(; x, y)
   T = EigenAnalysis(:VD)
   n, c = apply(T, t)
@@ -30,8 +30,8 @@
   @test Tables.matrix(t) ≈ Tables.matrix(tₒ)
 
   # SDS test
-  x = rand(Normal(0, 10), 1500)
-  y = x + rand(Normal(0, 2), 1500)
+  x = rand(rng, Normal(0, 10), 1500)
+  y = x + rand(rng, Normal(0, 2), 1500)
   t = Table(; x, y)
   T = EigenAnalysis(:VDV)
   n, c = apply(T, t)
@@ -73,8 +73,8 @@
   end
 
   # row table
-  x = rand(Normal(0, 10), 1500)
-  y = x + rand(Normal(0, 2), 1500)
+  x = rand(rng, Normal(0, 10), 1500)
+  y = x + rand(rng, Normal(0, 2), 1500)
   t = Table(; x, y)
   rt = Tables.rowtable(t)
   T = EigenAnalysis(:V)
@@ -84,9 +84,9 @@
   @test Tables.matrix(rt) ≈ Tables.matrix(rtₒ)
 
   # maxdim
-  x = randn(1000)
-  y = x + randn(1000)
-  z = 2x - y + randn(1000)
+  x = randn(rng, 1000)
+  y = x + randn(rng, 1000)
+  z = 2x - y + randn(rng, 1000)
   t = Table(; x, y, z)
 
   # PCA
diff --git a/test/transforms/filter.jl b/test/transforms/filter.jl
index 63ddd55a..9c56f187 100644
--- a/test/transforms/filter.jl
+++ b/test/transforms/filter.jl
@@ -104,7 +104,7 @@
   @test Tables.isrowtable(n)
 
   # performance tests
-  trng = MersenneTwister(2) # test rng
+  trng = StableRNG(2) # test rng
   x = rand(trng, 100_000)
   y = rand(trng, 100_000)
   c = CoDaArray((a=rand(trng, 100_000), b=rand(trng, 100_000), c=rand(trng, 100_000)))
diff --git a/test/transforms/projectionpursuit.jl b/test/transforms/projectionpursuit.jl
index de64b831..98439085 100644
--- a/test/transforms/projectionpursuit.jl
+++ b/test/transforms/projectionpursuit.jl
@@ -1,7 +1,7 @@
 @testset "ProjectionPursuit" begin
   @test TT.parameters(ProjectionPursuit()) == (tol=1.0e-6, deg=5, perc=0.9, n=100)
 
-  rng = MersenneTwister(42)
+  rng = StableRNG(42)
   N = 10_000
   a = [2randn(rng, N ÷ 2) .+ 6; randn(rng, N ÷ 2)]
   b = [3randn(rng, N ÷ 2); 2randn(rng, N ÷ 2)]
@@ -9,15 +9,15 @@
   d = c .+ 0.6randn(rng, N)
   t = (; a, b, c, d)
 
-  T = ProjectionPursuit(rng=MersenneTwister(2))
+  T = ProjectionPursuit(rng=StableRNG(2))
   n, c = apply(T, t)
 
   @test Tables.columnnames(n) == (:PP1, :PP2, :PP3, :PP4)
 
   μ = mean(Tables.matrix(n), dims=1)
   Σ = cov(Tables.matrix(n))
-  @test all(isapprox.(μ, 0, atol=1e-8))
-  @test all(isapprox.(Σ, I(4), atol=1e-8))
+  @test all(isapprox.(μ, 0, atol=1e-3))
+  @test all(isapprox.(Σ, I(4), atol=1e-2))
 
   tₒ = revert(T, n, c)
 
@@ -35,7 +35,7 @@
   b = rand(rng, BetaPrime(2), 4000)
   t = Table(; a, b)
 
-  T = ProjectionPursuit(rng=MersenneTwister(2))
+  T = ProjectionPursuit(rng=StableRNG(2))
   n, c = apply(T, t)
 
   μ = mean(Tables.matrix(n), dims=1)
diff --git a/test/transforms/sample.jl b/test/transforms/sample.jl
index 89a22991..5304ea2f 100644
--- a/test/transforms/sample.jl
+++ b/test/transforms/sample.jl
@@ -28,14 +28,14 @@
   @test n.c ⊆ t.c
   @test Tables.rowtable(n) == trows
 
-  T = Sample(8, replace=true, rng=MersenneTwister(2))
+  T = Sample(8, replace=true, rng=StableRNG(2))
   n, c = apply(T, t)
-  @test n.a == [3, 7, 8, 2, 2, 6, 2, 6]
-  @test n.b == [8, 2, 3, 1, 1, 5, 1, 5]
-  @test n.c == [1, 2, 9, 5, 5, 8, 5, 8]
+  @test n.a == [3, 3, 6, 3, 6, 6, 7, 3]
+  @test n.b == [4, 8, 5, 4, 5, 5, 2, 8]
+  @test n.c == [4, 1, 8, 4, 8, 8, 2, 1]
 
   w = pweights([0.1, 0.25, 0.15, 0.25, 0.1, 0.15])
-  T = Sample(10_000, w, replace=true, rng=MersenneTwister(2))
+  T = Sample(10_000, w, replace=true, rng=StableRNG(2))
   n, c = apply(T, t)
   nrows = Tables.rowtable(n)
   @test isapprox(count(==(trows[1]), nrows) / 10_000, 0.10, atol=0.01)
@@ -46,7 +46,7 @@
   @test isapprox(count(==(trows[6]), nrows) / 10_000, 0.15, atol=0.01)
 
   w = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0]
-  T = Sample(10_000, w, replace=true, rng=MersenneTwister(2))
+  T = Sample(10_000, w, replace=true, rng=StableRNG(2))
   n, c = apply(T, t)
   nrows = Tables.rowtable(n)
   @test isapprox(count(==(trows[1]), nrows) / 10_000, 1 / 21, atol=0.01)
@@ -57,7 +57,7 @@
   @test isapprox(count(==(trows[6]), nrows) / 10_000, 6 / 21, atol=0.01)
 
   # performance tests
-  trng = MersenneTwister(2) # test rng
+  trng = StableRNG(2) # test rng
   x = rand(trng, 100_000)
   y = rand(trng, 100_000)
   c = CoDaArray((a=rand(trng, 100_000), b=rand(trng, 100_000), c=rand(trng, 100_000)))