From f219195ab1acfb782648085bff58d169e7675322 Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Thu, 25 Jul 2024 00:34:58 -0400
Subject: [PATCH 01/16] Test scalar rand separately from vector rand

---
 test/fit.jl                                   |  2 +-
 test/multivariate/mvlognormal.jl              |  4 +-
 test/runtests.jl                              |  6 +-
 test/testutils.jl                             | 91 ++++++++++++++-----
 test/univariate/continuous.jl                 | 20 ----
 test/univariate/continuous/exponential.jl     | 11 +++
 test/univariate/continuous/inversegaussian.jl | 33 +++++++
 test/univariate/continuous/logitnormal.jl     |  9 ++
 test/univariate/continuous/lognormal.jl       | 14 +++
 test/univariate/continuous/normalcanon.jl     | 10 ++
 test/univariate/continuous/pareto.jl          | 10 ++
 .../continuous/pgeneralizedgaussian.jl        |  4 +-
 12 files changed, 165 insertions(+), 49 deletions(-)
 create mode 100644 test/univariate/continuous/inversegaussian.jl
 create mode 100644 test/univariate/continuous/normalcanon.jl
 create mode 100644 test/univariate/continuous/pareto.jl

diff --git a/test/fit.jl b/test/fit.jl
index 4483dd1ec6..b3474c6fe8 100644
--- a/test/fit.jl
+++ b/test/fit.jl
@@ -369,7 +369,7 @@ end
     for func in funcs, dist in (Laplace, Laplace{Float64})
         d = fit(dist, func[2](dist(5.0, 3.0), N + 1))
         @test isa(d, dist)
-        @test isapprox(location(d), 5.0, atol=0.02)
+        @test isapprox(location(d), 5.0, atol=0.021)
         @test isapprox(scale(d)   , 3.0, atol=0.03)
     end
 end
diff --git a/test/multivariate/mvlognormal.jl b/test/multivariate/mvlognormal.jl
index 89f923b889..b1c7f65758 100644
--- a/test/multivariate/mvlognormal.jl
+++ b/test/multivariate/mvlognormal.jl
@@ -105,8 +105,8 @@ end
     @test  entropy(l1)   ≈ entropy(l2)
     @test logpdf(l1,5.0) ≈ logpdf(l2,[5.0])
     @test pdf(l1,5.0)    ≈ pdf(l2,[5.0])
-    @test (Random.seed!(78393) ; [rand(l1)]) == (Random.seed!(78393) ; rand(l2))
-    @test [rand(MersenneTwister(78393), l1)] == rand(MersenneTwister(78393), l2)
+    @test (Random.seed!(78393) ; [rand(l1)]) ≈ (Random.seed!(78393) ; rand(l2))
+    @test [rand(MersenneTwister(78393), l1)] ≈ rand(MersenneTwister(78393), l2)
 end
 
 ###### General Testing
diff --git a/test/runtests.jl b/test/runtests.jl
index 583132c536..48758e82f0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -22,6 +22,7 @@ const tests = [
     "truncated/discrete_uniform",
     "censored",
     "univariate/continuous/normal",
+    "univariate/continuous/normalcanon",
     "univariate/continuous/laplace",
     "univariate/continuous/cauchy",
     "univariate/continuous/uniform",
@@ -77,12 +78,14 @@ const tests = [
     "univariate/continuous/exponential",
     "univariate/continuous/gamma",
     "univariate/continuous/gumbel",
+    "univariate/continuous/inversegaussian",
     "univariate/continuous/lindley",
     "univariate/continuous/logistic",
     "univariate/continuous/johnsonsu",
     "univariate/continuous/noncentralchisq",
     "univariate/continuous/weibull",
     "pdfnorm",
+    "univariate/continuous/pareto",
     "univariate/continuous/rician",
     "functionals",
     "density_interface",
@@ -137,15 +140,12 @@ const tests = [
     # "univariate/continuous/generalizedextremevalue",
     # "univariate/continuous/generalizedpareto",
     # "univariate/continuous/inversegamma",
-    # "univariate/continuous/inversegaussian",
     # "univariate/continuous/ksdist",
     # "univariate/continuous/ksonesided",
     # "univariate/continuous/levy",
     # "univariate/continuous/noncentralbeta",
     # "univariate/continuous/noncentralf",
-    # "univariate/continuous/normalcanon",
     # "univariate/continuous/normalinversegaussian",
-    # "univariate/continuous/pareto",
     # "univariate/continuous/rayleigh",
     # "univariate/continuous/studentizedrange",
     # "univariate/continuous/symtriangular",
diff --git a/test/testutils.jl b/test/testutils.jl
index 2859856a73..358eb6c952 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -28,8 +28,9 @@ end
 # testing the implementation of a discrete univariate distribution
 #
 function test_distr(distr::DiscreteUnivariateDistribution, n::Int;
-                    testquan::Bool=true)
-
+                    testquan::Bool=true, rng::AbstractRNG=MersenneTwister(),
+                    test_scalar_rand::Bool=false)
+    println("    testing $(distr)")
     test_range(distr)
     vs = get_evalsamples(distr, 0.00001)
 
@@ -40,7 +41,12 @@ function test_distr(distr::DiscreteUnivariateDistribution, n::Int;
 
     test_stats(distr, vs)
     test_samples(distr, n)
-    test_samples(distr, n, rng=MersenneTwister())
+    test_samples(distr, n, rng=rng)
+    if test_scalar_rand
+        xs = test_samples(distr, n; call_scalar = true)
+        xs = test_samples(distr, n, rng=rng, call_scalar = true)
+    end
+
     test_params(distr)
 end
 
@@ -67,7 +73,9 @@ end
 # testing the implementation of a continuous univariate distribution
 #
 function test_distr(distr::ContinuousUnivariateDistribution, n::Int;
-                    testquan::Bool=true, rng::AbstractRNG=MersenneTwister(123))
+                    testquan::Bool=true, rng::AbstractRNG=MersenneTwister(123),
+                    test_scalar_rand::Bool=false)
+    println("    testing $(distr)")
     test_range(distr)
     vs = get_evalsamples(distr, 0.01, 2000)
 
@@ -82,6 +90,13 @@ function test_distr(distr::ContinuousUnivariateDistribution, n::Int;
     allow_test_stats(distr) && test_stats(distr, xs)
     xs = test_samples(distr, n, rng=rng)
     allow_test_stats(distr) && test_stats(distr, xs)
+    if test_scalar_rand
+        xs = test_samples(distr, n; call_scalar = true)
+        allow_test_stats(distr) && test_stats(distr, xs)
+        xs = test_samples(distr, n, rng=rng, call_scalar = true)
+        allow_test_stats(distr) && test_stats(distr, xs)
+    end
+
     test_params(distr)
 end
 
@@ -101,8 +116,8 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
                       n::Int;                                   # number of samples to generate
                       q::Float64=1.0e-7,                        # confidence interval, 1 - q as confidence
                       verbose::Bool=false,                      # show intermediate info (for debugging)
-                      rng::Union{AbstractRNG, Missing}=missing) # add an rng?
-
+                      rng::Union{AbstractRNG, Missing}=missing,  # add an rng?
+                      call_scalar::Bool=false)                  # directly scall the scalar rand(d) instead
     # The basic idea
     # ------------------
     #   Generate n samples, and count the occurrences of each value within a reasonable range.
@@ -146,14 +161,30 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
     # generate samples using RNG passed or default RNG
     # we also check reproducibility
     if rng === missing
+        Random.seed!(1234)    
+        samples = if !call_scalar
+            rand(s, n)
+        else
+            map((_) -> rand(s), 1:n)
+        end
         Random.seed!(1234)
-        samples = rand(s, n)
-        Random.seed!(1234)
-        samples2 = rand(s, n)
+        samples2 = if !call_scalar
+            rand(s, n)
+        else
+            map((_) -> rand(s), 1:n)
+        end
     else
         rng2 = deepcopy(rng)
-        samples = rand(rng, s, n)
-        samples2 = rand(rng2, s, n)
+        samples = if !call_scalar
+            rand(rng, s, n)
+        else
+            map((_) -> rand(rng, s), 1:n)
+        end
+        samples2 = if !call_scalar
+            rand(rng2, s, n)
+        else
+            map((_) -> rand(rng2, s), 1:n)
+        end
     end
     @test length(samples) == n
     @test samples2 == samples
@@ -180,8 +211,8 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
 end
 
 test_samples(distr::DiscreteUnivariateDistribution, n::Int;
-             q::Float64=1.0e-6, verbose::Bool=false, rng=missing) =
-    test_samples(distr, distr, n; q=q, verbose=verbose, rng=rng)
+             q::Float64=1.0e-6, verbose::Bool=false, rng=missing, kwargs...) =
+    test_samples(distr, distr, n; q=q, verbose=verbose, rng=rng, kwargs...)
 
 # for continuous samplers
 #
@@ -191,7 +222,8 @@ function test_samples(s::Sampleable{Univariate, Continuous},    # the sampleable
                       nbins::Int=50,                            # divide the main interval into nbins
                       q::Float64=1.0e-6,                        # confidence interval, 1 - q as confidence
                       verbose::Bool=false,                      # show intermediate info (for debugging)
-                      rng::Union{AbstractRNG, Missing}=missing) # add an rng?
+                      rng::Union{AbstractRNG, Missing}=missing, # add an rng?
+                      call_scalar::Bool=false)                  # directly scall the scalar rand(d) instead
 
     # The basic idea
     # ------------------
@@ -246,14 +278,30 @@ function test_samples(s::Sampleable{Univariate, Continuous},    # the sampleable
     # generate samples using RNG passed or default RNG
     # we also check reproducibility
     if rng === missing
+        Random.seed!(1234)    
+        samples = if !call_scalar
+            rand(s, n)
+        else
+            map((_) -> rand(s), 1:n)
+        end
         Random.seed!(1234)
-        samples = rand(s, n)
-        Random.seed!(1234)
-        samples2 = rand(s, n)
+        samples2 = if !call_scalar
+            rand(s, n)
+        else
+            map((_) -> rand(s), 1:n)
+        end
     else
         rng2 = deepcopy(rng)
-        samples = rand(rng, s, n)
-        samples2 = rand(rng2, s, n)
+        samples = if !call_scalar
+            rand(rng, s, n)
+        else
+            map((_) -> rand(rng, s), 1:n)
+        end
+        samples2 = if !call_scalar
+            rand(rng2, s, n)
+        else
+            map((_) -> rand(rng2, s), 1:n)
+        end
     end
     @test length(samples) == n
     @test samples2 == samples
@@ -284,8 +332,8 @@ function test_samples(s::Sampleable{Univariate, Continuous},    # the sampleable
     return samples
 end
 
-test_samples(distr::ContinuousUnivariateDistribution, n::Int; nbins::Int=50, q::Float64=1.0e-6, verbose::Bool=false, rng=missing) =
-    test_samples(distr, distr, n; nbins=nbins, q=q, verbose=verbose, rng=rng)
+test_samples(distr::ContinuousUnivariateDistribution, n::Int; nbins::Int=50, q::Float64=1.0e-6, verbose::Bool=false, rng=missing, kwargs...) =
+    test_samples(distr, distr, n; nbins=nbins, q=q, verbose=verbose, rng=rng, kwargs...)
 
 
 #### Testing range & support methods
@@ -583,6 +631,7 @@ end
 allow_test_stats(d::UnivariateDistribution) = true
 allow_test_stats(d::NoncentralBeta) = false
 allow_test_stats(::StudentizedRange) = false
+allow_test_stats(::LogitNormal) = false
 
 function test_stats(d::ContinuousUnivariateDistribution, xs::AbstractVector{Float64})
     # using Monte Carlo methods
diff --git a/test/univariate/continuous.jl b/test/univariate/continuous.jl
index 3a1da94b92..fcaf91b7a6 100644
--- a/test/univariate/continuous.jl
+++ b/test/univariate/continuous.jl
@@ -79,26 +79,6 @@ end
     @test cdf(d, 1.0) ≈ 0.000787319 atol=1e-9
 end
 
-@testset "NormalInverseGaussian random repeatable and basic metrics" begin
-    rng = Random.MersenneTwister(42)
-    rng2 = copy(rng)
-    µ = 0.0
-    α = 1.0
-    β = 0.5
-    δ = 3.0
-    g = sqrt(α^2 - β^2)
-    d = NormalInverseGaussian(μ, α, β, δ)
-    v1 = rand(rng, d)
-    v2 = rand(rng, d)
-    v3 = rand(rng2, d)
-    @test v1 ≈ v3
-    @test v1 ≉ v2
-
-    @test mean(d) ≈ µ + β * δ / g
-    @test var(d) ≈ δ * α^2 / g^3
-    @test skewness(d) ≈ 3β/(α*sqrt(δ*g))
-end
-
 @testset "edge cases" begin
     # issue #1371: cdf should not return -0.0
     @test cdf(Rayleigh(1), 0) === 0.0
diff --git a/test/univariate/continuous/exponential.jl b/test/univariate/continuous/exponential.jl
index 6704554a79..90e4787384 100644
--- a/test/univariate/continuous/exponential.jl
+++ b/test/univariate/continuous/exponential.jl
@@ -2,3 +2,14 @@
 test_cgf(Exponential(1), (0.9, -1, -100f0, -1e6))
 test_cgf(Exponential(0.91), (0.9, -1, -100f0, -1e6))
 test_cgf(Exponential(10  ), (0.08, -1, -100f0, -1e6))
+
+# Sampling Tests
+@testset "Exponential sampling tests" begin
+    for d in [
+        Exponential(1),
+        Exponential(0.91),
+        Exponential(10)
+    ]
+        test_distr(d, 10^6, test_scalar_rand = true)
+    end
+end
\ No newline at end of file
diff --git a/test/univariate/continuous/inversegaussian.jl b/test/univariate/continuous/inversegaussian.jl
new file mode 100644
index 0000000000..1c28202af8
--- /dev/null
+++ b/test/univariate/continuous/inversegaussian.jl
@@ -0,0 +1,33 @@
+@testset "NormalInverseGaussian random repeatable and basic metrics" begin
+    rng = Random.MersenneTwister(42)
+    rng2 = copy(rng)
+    µ = 0.0
+    α = 1.0
+    β = 0.5
+    δ = 3.0
+    g = sqrt(α^2 - β^2)
+    d = NormalInverseGaussian(μ, α, β, δ)
+    v1 = rand(rng, d)
+    v2 = rand(rng, d)
+    v3 = rand(rng2, d)
+    @test v1 ≈ v3
+    @test v1 ≉ v2
+
+    @test mean(d) ≈ µ + β * δ / g
+    @test var(d) ≈ δ * α^2 / g^3
+    @test skewness(d) ≈ 3β/(α*sqrt(δ*g))
+end
+
+# Sampling Tests
+@testset "InverseGaussian sampling tests" begin
+    for d in [
+        InverseGaussian()
+        InverseGaussian(0.8)
+        InverseGaussian(2.0)
+        InverseGaussian(1.0, 1.0)
+        InverseGaussian(2.0, 1.5)
+        InverseGaussian(2.0, 7.0)
+    ]
+        test_distr(d, 10^6, test_scalar_rand = true)
+    end
+end
\ No newline at end of file
diff --git a/test/univariate/continuous/logitnormal.jl b/test/univariate/continuous/logitnormal.jl
index 454376606c..b8c4a578d8 100644
--- a/test/univariate/continuous/logitnormal.jl
+++ b/test/univariate/continuous/logitnormal.jl
@@ -64,3 +64,12 @@ end
     @test convert(LogitNormal{Float32}, d) === d
     @test typeof(convert(LogitNormal{Float64}, d)) == typeof(LogitNormal(2,1))
 end
+
+@testset "Logitnormal Sampling Tests" begin
+    for d in [
+        LogitNormal(-2, 3),
+        LogitNormal(0, 0.2)
+    ]
+        test_distr(d, 10^6, test_scalar_rand = true)
+    end
+end
\ No newline at end of file
diff --git a/test/univariate/continuous/lognormal.jl b/test/univariate/continuous/lognormal.jl
index 4a27f2e8ac..26c16f7730 100644
--- a/test/univariate/continuous/lognormal.jl
+++ b/test/univariate/continuous/lognormal.jl
@@ -314,3 +314,17 @@ end
     @test @inferred(gradlogpdf(LogNormal(0.0, 1.0), BigFloat(-1))) == big(0.0)
     @test isnan_type(BigFloat, @inferred(gradlogpdf(LogNormal(0.0, 1.0), BigFloat(NaN))))
 end
+
+@testset "LogNormal Sampling Tests" begin
+    for d in [
+            LogNormal()
+            LogNormal(1.0)
+            LogNormal(0.0, 2.0)
+            LogNormal(1.0, 2.0)
+            LogNormal(3.0, 0.5)
+            LogNormal(3.0, 1.0)
+            LogNormal(3.0, 2.0)
+        ]
+        test_distr(d, 10^6, test_scalar_rand = true)
+    end
+end
\ No newline at end of file
diff --git a/test/univariate/continuous/normalcanon.jl b/test/univariate/continuous/normalcanon.jl
new file mode 100644
index 0000000000..443d6ec494
--- /dev/null
+++ b/test/univariate/continuous/normalcanon.jl
@@ -0,0 +1,10 @@
+# Sampling Tests
+@testset "NormalCanon sampling tests" begin
+    for d in [
+        NormalCanon()
+        NormalCanon(-1.0, 2.5)
+        NormalCanon(2.0, 0.8)
+    ]
+        test_distr(d, 10^6, test_scalar_rand = true)
+    end
+end
\ No newline at end of file
diff --git a/test/univariate/continuous/pareto.jl b/test/univariate/continuous/pareto.jl
new file mode 100644
index 0000000000..9e7e95775e
--- /dev/null
+++ b/test/univariate/continuous/pareto.jl
@@ -0,0 +1,10 @@
+@testset "Pareto Sampling Tests" begin
+    for d in [
+            Pareto()
+            Pareto(2.0)
+            Pareto(2.0, 1.5)
+            Pareto(3.0, 2.0)
+        ]
+        test_distr(d, 10^6, test_scalar_rand = true)
+    end
+end
\ No newline at end of file
diff --git a/test/univariate/continuous/pgeneralizedgaussian.jl b/test/univariate/continuous/pgeneralizedgaussian.jl
index 49de167cc5..52ba0ee6b0 100644
--- a/test/univariate/continuous/pgeneralizedgaussian.jl
+++ b/test/univariate/continuous/pgeneralizedgaussian.jl
@@ -74,7 +74,7 @@ using Test
             end
 
             # Additional tests, including sampling
-            test_distr(d, 10^6)
+            test_distr(d, 10^6, test_scalar_rand = true)
         end
     end
 
@@ -104,6 +104,6 @@ using Test
         @test quantile(d, 1 // 2) ≈ μ
 
         # Additional tests, including sampling
-        test_distr(d, 10^6)
+        test_distr(d, 10^6, test_scalar_rand = true)
     end
 end

From 01cba644ee71261c776d065123081ac6a0f9c25c Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Thu, 25 Jul 2024 03:19:48 -0400
Subject: [PATCH 02/16] Add specialized rand! for many distributions

---
 src/univariate/continuous/exponential.jl        |  2 ++
 src/univariate/continuous/inversegaussian.jl    | 17 ++++++++++++++++-
 src/univariate/continuous/laplace.jl            |  8 ++++++++
 src/univariate/continuous/logitnormal.jl        |  6 +++++-
 src/univariate/continuous/lognormal.jl          |  5 ++++-
 src/univariate/continuous/normal.jl             |  2 +-
 src/univariate/continuous/normalcanon.jl        |  6 ++++++
 src/univariate/continuous/pareto.jl             |  6 ++++++
 .../continuous/pgeneralizedgaussian.jl          | 13 ++++++++++++-
 9 files changed, 60 insertions(+), 5 deletions(-)

diff --git a/src/univariate/continuous/exponential.jl b/src/univariate/continuous/exponential.jl
index d1c487057f..a14394d564 100644
--- a/src/univariate/continuous/exponential.jl
+++ b/src/univariate/continuous/exponential.jl
@@ -107,6 +107,8 @@ cf(d::Exponential, t::Real) = 1/(1 - t * im * scale(d))
 #### Sampling
 rand(rng::AbstractRNG, d::Exponential) = xval(d, randexp(rng))
 
+rand!(rng::AbstractRNG, d::Exponential, A::AbstractArray{<:Real}) =
+    A .= xval.(d, randexp!(rng, A))
 
 #### Fit model
 
diff --git a/src/univariate/continuous/inversegaussian.jl b/src/univariate/continuous/inversegaussian.jl
index a076d85b21..e48b28bc74 100644
--- a/src/univariate/continuous/inversegaussian.jl
+++ b/src/univariate/continuous/inversegaussian.jl
@@ -163,8 +163,23 @@ function rand(rng::AbstractRNG, d::InverseGaussian)
     u >= p1 ? μ^2 / x1 : x1
 end
 
-#### Fit model
+function rand!(rng::AbstractRNG, d::InverseGaussian, A::AbstractArray{<:Real})
+    # Based off of the non-vectorized code
+    μ, λ = params(d)
+    Z = V = W = X1 = A # prevents extra heap allocs
+    randn!(rng, A)
+    V .*= Z
+    W .*= μ
+    X1 .= μ .+ μ / (2λ) .* (W .- sqrt.(W .* (4λ .+ W)))
+    A .= ifelse.( # TODO Can we avoid heap-allocing a whole float vec?
+        rand(rng, size(A)...) .>= μ ./ (μ .+ X1),
+        (μ * μ) ./ X1,
+        X1
+    )
+end
 
+
+#### Fit model
 """
 Sufficient statistics for `InverseGaussian`, containing the weighted
 sum of observations, the weighted sum of inverse points and sum of weights.
diff --git a/src/univariate/continuous/laplace.jl b/src/univariate/continuous/laplace.jl
index 2a7bf04a47..5bb6417ad7 100644
--- a/src/univariate/continuous/laplace.jl
+++ b/src/univariate/continuous/laplace.jl
@@ -126,6 +126,14 @@ Base.:*(c::Real, d::Laplace) = Laplace(c * d.μ, abs(c) * d.θ)
 rand(rng::AbstractRNG, d::Laplace) =
     d.μ + d.θ*randexp(rng)*ifelse(rand(rng, Bool), 1, -1)
 
+function rand!(rng::AbstractRNG, d::Laplace, A::AbstractArray{<:Real})
+    randexp!(rng, A)
+    A .= muladd.(
+        A,
+        ifelse.(rand(rng, Bool, size(A)), d.θ, -d.θ),
+        d.μ
+    )
+end
 
 #### Fitting
 
diff --git a/src/univariate/continuous/logitnormal.jl b/src/univariate/continuous/logitnormal.jl
index 7fd5fb9113..5d26f93721 100644
--- a/src/univariate/continuous/logitnormal.jl
+++ b/src/univariate/continuous/logitnormal.jl
@@ -157,7 +157,11 @@ end
 
 #### Sampling
 
-rand(rng::AbstractRNG, d::LogitNormal) = logistic(randn(rng) * d.σ + d.μ)
+rand(rng::AbstractRNG, d::LogitNormal) = logistic(muladd(d.σ, randn(rng), d.μ))
+
+rand!(rng::AbstractRNG, d::LogitNormal, A::AbstractArray{<:Real}) =
+    A .= logistic.(muladd.(d.σ, randn!(rng, A), d.μ))
+
 
 ## Fitting
 
diff --git a/src/univariate/continuous/lognormal.jl b/src/univariate/continuous/lognormal.jl
index 666b14248a..5c5e529884 100644
--- a/src/univariate/continuous/lognormal.jl
+++ b/src/univariate/continuous/lognormal.jl
@@ -156,7 +156,10 @@ end
 
 #### Sampling
 
-rand(rng::AbstractRNG, d::LogNormal) = exp(randn(rng) * d.σ + d.μ)
+rand(rng::AbstractRNG, d::LogNormal) = exp(muladd(d.σ, randn(rng), d.μ))
+
+rand!(rng::AbstractRNG, d::LogNormal, A::AbstractArray{<:Real}) =
+    A .= exp.(muladd.(d.σ, randn!(rng, A), d.μ))
 
 ## Fitting
 
diff --git a/src/univariate/continuous/normal.jl b/src/univariate/continuous/normal.jl
index 4231522674..3e7d446b3e 100644
--- a/src/univariate/continuous/normal.jl
+++ b/src/univariate/continuous/normal.jl
@@ -114,7 +114,7 @@ Base.:*(c::Real, d::Normal) = Normal(c * d.μ, abs(c) * d.σ)
 
 #### Sampling
 
-rand(rng::AbstractRNG, d::Normal{T}) where {T} = d.μ + d.σ * randn(rng, float(T))
+rand(rng::AbstractRNG, d::Normal{T}) where {T} = muladd(d.σ, randn(rng, float(T)), d.μ)
 
 rand!(rng::AbstractRNG, d::Normal, A::AbstractArray{<:Real}) = A .= muladd.(d.σ, randn!(rng, A), d.μ)
 
diff --git a/src/univariate/continuous/normalcanon.jl b/src/univariate/continuous/normalcanon.jl
index de29c4b196..0272f1c05c 100644
--- a/src/univariate/continuous/normalcanon.jl
+++ b/src/univariate/continuous/normalcanon.jl
@@ -89,6 +89,12 @@ invlogccdf(d::NormalCanon, lp::Real) = xval(d, norminvlogccdf(lp))
 
 rand(rng::AbstractRNG, cf::NormalCanon) = cf.μ + randn(rng) / sqrt(cf.λ)
 
+rand!(rng::AbstractRNG, cf::NormalCanon, A::AbstractArray{<:Real}) =
+    # A .* inv(sqrt(cf.λ)) is faster but less accurate than A ./ sqrt(cf.λ)
+    # in floating point. This is a fast-math style optimization.
+    # It shouldn't be a problem for PRNG.
+    A .= muladd.(randn!(rng, A), inv(sqrt(cf.λ)), cf.μ)
+
 #### Affine transformations
 
 function Base.:+(d::NormalCanon, c::Real)
diff --git a/src/univariate/continuous/pareto.jl b/src/univariate/continuous/pareto.jl
index 10bb246c5f..0d2c0d024c 100644
--- a/src/univariate/continuous/pareto.jl
+++ b/src/univariate/continuous/pareto.jl
@@ -112,6 +112,12 @@ quantile(d::Pareto, p::Real) = cquantile(d, 1 - p)
 
 rand(rng::AbstractRNG, d::Pareto) = d.θ * exp(randexp(rng) / d.α)
 
+rand!(rng::AbstractRNG, d::Pareto, A::AbstractArray{<:Real}) =
+    # A .* 1/d.α is faster than A ./ d.α, but less exact
+    # fast-math style optimization; should be fine for PRNG
+    A .=  d.θ .* exp.(randexp!(rng, A) .* (1/d.α))
+
+
 ## Fitting
 
 function fit_mle(::Type{<:Pareto}, x::AbstractArray{T}) where T<:Real
diff --git a/src/univariate/continuous/pgeneralizedgaussian.jl b/src/univariate/continuous/pgeneralizedgaussian.jl
index 9d41d9d464..4c9dcfa5ff 100644
--- a/src/univariate/continuous/pgeneralizedgaussian.jl
+++ b/src/univariate/continuous/pgeneralizedgaussian.jl
@@ -141,9 +141,20 @@ function rand(rng::AbstractRNG, d::PGeneralizedGaussian)
     inv_p = inv(d.p)
     g = Gamma(inv_p, 1)
     z = d.α * rand(rng, g)^inv_p
-    if rand(rng) < 0.5
+    if rand(rng, Bool)
         return d.μ - z
     else
         return d.μ + z
     end
 end
+
+function rand!(rng::AbstractRNG, d::PGeneralizedGaussian, A::AbstractArray{<:Real})
+    inv_p = inv(d.p)
+    g = Gamma(inv_p, 1)
+    A .= rand!(rng, g, A).^inv_p
+    A .= muladd.(
+        A,
+        ifelse.(rand(rng, Bool, size(A)), d.α, -d.α),
+        d.μ
+    )
+end
\ No newline at end of file

From 7cf954c3dd7a5507480a02f69b9f648eda1fd9c4 Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Fri, 9 Aug 2024 03:46:41 -0400
Subject: [PATCH 03/16] Restore location of old NormalInverseGaussian tests

---
 test/univariate/continuous.jl                 | 20 +++++++++++++++++++
 test/univariate/continuous/inversegaussian.jl | 20 -------------------
 2 files changed, 20 insertions(+), 20 deletions(-)

diff --git a/test/univariate/continuous.jl b/test/univariate/continuous.jl
index fcaf91b7a6..3a1da94b92 100644
--- a/test/univariate/continuous.jl
+++ b/test/univariate/continuous.jl
@@ -79,6 +79,26 @@ end
     @test cdf(d, 1.0) ≈ 0.000787319 atol=1e-9
 end
 
+@testset "NormalInverseGaussian random repeatable and basic metrics" begin
+    rng = Random.MersenneTwister(42)
+    rng2 = copy(rng)
+    µ = 0.0
+    α = 1.0
+    β = 0.5
+    δ = 3.0
+    g = sqrt(α^2 - β^2)
+    d = NormalInverseGaussian(μ, α, β, δ)
+    v1 = rand(rng, d)
+    v2 = rand(rng, d)
+    v3 = rand(rng2, d)
+    @test v1 ≈ v3
+    @test v1 ≉ v2
+
+    @test mean(d) ≈ µ + β * δ / g
+    @test var(d) ≈ δ * α^2 / g^3
+    @test skewness(d) ≈ 3β/(α*sqrt(δ*g))
+end
+
 @testset "edge cases" begin
     # issue #1371: cdf should not return -0.0
     @test cdf(Rayleigh(1), 0) === 0.0
diff --git a/test/univariate/continuous/inversegaussian.jl b/test/univariate/continuous/inversegaussian.jl
index 1c28202af8..ecbb48d139 100644
--- a/test/univariate/continuous/inversegaussian.jl
+++ b/test/univariate/continuous/inversegaussian.jl
@@ -1,23 +1,3 @@
-@testset "NormalInverseGaussian random repeatable and basic metrics" begin
-    rng = Random.MersenneTwister(42)
-    rng2 = copy(rng)
-    µ = 0.0
-    α = 1.0
-    β = 0.5
-    δ = 3.0
-    g = sqrt(α^2 - β^2)
-    d = NormalInverseGaussian(μ, α, β, δ)
-    v1 = rand(rng, d)
-    v2 = rand(rng, d)
-    v3 = rand(rng2, d)
-    @test v1 ≈ v3
-    @test v1 ≉ v2
-
-    @test mean(d) ≈ µ + β * δ / g
-    @test var(d) ≈ δ * α^2 / g^3
-    @test skewness(d) ≈ 3β/(α*sqrt(δ*g))
-end
-
 # Sampling Tests
 @testset "InverseGaussian sampling tests" begin
     for d in [

From 3a607943c4c60eecf70a71fbdcc785820ac0ad3a Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Fri, 9 Aug 2024 14:35:19 -0400
Subject: [PATCH 04/16] Remove duplication of inversegaussian in runtests.jl

---
 test/runtests.jl | 1 -
 1 file changed, 1 deletion(-)

diff --git a/test/runtests.jl b/test/runtests.jl
index ee2952c3cd..cd3afe17f9 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -78,7 +78,6 @@ const tests = [
     "univariate/continuous/exponential",
     "univariate/continuous/gamma",
     "univariate/continuous/gumbel",
-    "univariate/continuous/inversegaussian",
     "univariate/continuous/lindley",
     "univariate/continuous/logistic",
     "univariate/continuous/johnsonsu",

From 368dc9c8748b8847087252869911641a20af8a60 Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Sat, 17 Aug 2024 20:03:30 -0400
Subject: [PATCH 05/16] Apply many suggestions from code review

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
---
 src/univariate/continuous/exponential.jl |  7 +++++--
 src/univariate/continuous/logitnormal.jl | 11 +++++++----
 src/univariate/continuous/lognormal.jl   | 10 +++++++---
 src/univariate/continuous/normal.jl      |  9 +++++++--
 src/univariate/continuous/normalcanon.jl | 10 +++++-----
 src/univariate/continuous/pareto.jl      | 13 +++++++------
 6 files changed, 38 insertions(+), 22 deletions(-)

diff --git a/src/univariate/continuous/exponential.jl b/src/univariate/continuous/exponential.jl
index a14394d564..53d8161d0d 100644
--- a/src/univariate/continuous/exponential.jl
+++ b/src/univariate/continuous/exponential.jl
@@ -107,8 +107,11 @@ cf(d::Exponential, t::Real) = 1/(1 - t * im * scale(d))
 #### Sampling
 rand(rng::AbstractRNG, d::Exponential) = xval(d, randexp(rng))
 
-rand!(rng::AbstractRNG, d::Exponential, A::AbstractArray{<:Real}) =
-    A .= xval.(d, randexp!(rng, A))
+function rand!(rng::AbstractRNG, d::Exponential, A::AbstractArray{<:Real})
+    randexp!(rng, A)
+    map!(Base.Fix1(xval, d), A, A)
+    return A
+end
 
 #### Fit model
 
diff --git a/src/univariate/continuous/logitnormal.jl b/src/univariate/continuous/logitnormal.jl
index 5d26f93721..cf0832d00c 100644
--- a/src/univariate/continuous/logitnormal.jl
+++ b/src/univariate/continuous/logitnormal.jl
@@ -157,11 +157,14 @@ end
 
 #### Sampling
 
-rand(rng::AbstractRNG, d::LogitNormal) = logistic(muladd(d.σ, randn(rng), d.μ))
-
-rand!(rng::AbstractRNG, d::LogitNormal, A::AbstractArray{<:Real}) =
-    A .= logistic.(muladd.(d.σ, randn!(rng, A), d.μ))
+xval(d::LogitNormal, z::Real) = logistic(muladd(d.σ, z, d.μ))
 
+rand(rng::AbstractRNG, d::LogitNormal) = xval(d, randn(rng))
+function rand!(rng::AbstractRNG, d::LogitNormal, A::AbstractArray{<:Real})
+    randn!(rng, A)
+    map!(Base.Fix1(xval, d), A, A)
+    return A
+end
 
 ## Fitting
 
diff --git a/src/univariate/continuous/lognormal.jl b/src/univariate/continuous/lognormal.jl
index 5c5e529884..6969b4c9d8 100644
--- a/src/univariate/continuous/lognormal.jl
+++ b/src/univariate/continuous/lognormal.jl
@@ -156,10 +156,14 @@ end
 
 #### Sampling
 
-rand(rng::AbstractRNG, d::LogNormal) = exp(muladd(d.σ, randn(rng), d.μ))
+xval(d::LogNormal, z::Real) = exp(muladd(d.σ, z, d.μ)  
 
-rand!(rng::AbstractRNG, d::LogNormal, A::AbstractArray{<:Real}) =
-    A .= exp.(muladd.(d.σ, randn!(rng, A), d.μ))
+rand(rng::AbstractRNG, d::LogNormal) = xval(d,  randn(rng))
+function rand!(rng::AbstractRNG, d::LogNormal, A::AbstractArray{<:Real})
+    randn!(rng, A)
+    map!(Base.Fix1(xval, d), A, A)
+    return A
+end
 
 ## Fitting
 
diff --git a/src/univariate/continuous/normal.jl b/src/univariate/continuous/normal.jl
index 3e7d446b3e..ef2b77bb58 100644
--- a/src/univariate/continuous/normal.jl
+++ b/src/univariate/continuous/normal.jl
@@ -114,9 +114,14 @@ Base.:*(c::Real, d::Normal) = Normal(c * d.μ, abs(c) * d.σ)
 
 #### Sampling
 
-rand(rng::AbstractRNG, d::Normal{T}) where {T} = muladd(d.σ, randn(rng, float(T)), d.μ)
+xval(d::Normal, z::Real) = muladd(d.σ, z, d.μ)
 
-rand!(rng::AbstractRNG, d::Normal, A::AbstractArray{<:Real}) = A .= muladd.(d.σ, randn!(rng, A), d.μ)
+rand(rng::AbstractRNG, d::Normal{T}) where {T} = xval(d, randn(rng, float(T)))
+function rand!(rng::AbstractRNG, d::Normal, A::AbstractArray{<:Real})
+    randn!(rng, A)
+    map!(Base.Fix1(xval, d), A, A)
+    return A
+end
 
 #### Fitting
 
diff --git a/src/univariate/continuous/normalcanon.jl b/src/univariate/continuous/normalcanon.jl
index 0272f1c05c..36f1215be2 100644
--- a/src/univariate/continuous/normalcanon.jl
+++ b/src/univariate/continuous/normalcanon.jl
@@ -89,11 +89,11 @@ invlogccdf(d::NormalCanon, lp::Real) = xval(d, norminvlogccdf(lp))
 
 rand(rng::AbstractRNG, cf::NormalCanon) = cf.μ + randn(rng) / sqrt(cf.λ)
 
-rand!(rng::AbstractRNG, cf::NormalCanon, A::AbstractArray{<:Real}) =
-    # A .* inv(sqrt(cf.λ)) is faster but less accurate than A ./ sqrt(cf.λ)
-    # in floating point. This is a fast-math style optimization.
-    # It shouldn't be a problem for PRNG.
-    A .= muladd.(randn!(rng, A), inv(sqrt(cf.λ)), cf.μ)
+function rand!(rng::AbstractRNG, cf::NormalCanon, A::AbstractArray{<:Real})
+    randn!(rng, A)
+    A .= cf.μ .+ A ./ sqrt(cf.λ)
+    return A
+end
 
 #### Affine transformations
 
diff --git a/src/univariate/continuous/pareto.jl b/src/univariate/continuous/pareto.jl
index 0d2c0d024c..249fa29eb7 100644
--- a/src/univariate/continuous/pareto.jl
+++ b/src/univariate/continuous/pareto.jl
@@ -110,13 +110,14 @@ quantile(d::Pareto, p::Real) = cquantile(d, 1 - p)
 
 #### Sampling
 
-rand(rng::AbstractRNG, d::Pareto) = d.θ * exp(randexp(rng) / d.α)
-
-rand!(rng::AbstractRNG, d::Pareto, A::AbstractArray{<:Real}) =
-    # A .* 1/d.α is faster than A ./ d.α, but less exact
-    # fast-math style optimization; should be fine for PRNG
-    A .=  d.θ .* exp.(randexp!(rng, A) .* (1/d.α))
+xval(d::Pareto, z::Real) = d.θ * exp(z / d.α)
 
+rand(rng::AbstractRNG, d::Pareto) = xval(d, randexp(rng))
+function rand!(rng::AbstractRNG, d::Pareto, A::AbstractArray{<:Real})
+    randexp!(rng, A)
+    map!(Base.Fix1(xval, d), A, A)
+    return A
+end
 
 ## Fitting
 

From ad61b583c59cc45cbebf2096ad18c6c7141fcbfc Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Sat, 17 Aug 2024 20:14:00 -0400
Subject: [PATCH 06/16] Apply other suggestions

---
 src/univariate/continuous/inversegaussian.jl     | 16 ----------------
 src/univariate/continuous/laplace.jl             |  9 ---------
 src/univariate/continuous/lognormal.jl           |  2 +-
 .../continuous/pgeneralizedgaussian.jl           | 11 -----------
 4 files changed, 1 insertion(+), 37 deletions(-)

diff --git a/src/univariate/continuous/inversegaussian.jl b/src/univariate/continuous/inversegaussian.jl
index 2cda7730ca..ec507bc610 100644
--- a/src/univariate/continuous/inversegaussian.jl
+++ b/src/univariate/continuous/inversegaussian.jl
@@ -169,22 +169,6 @@ function rand(rng::AbstractRNG, d::InverseGaussian)
     u >= p1 ? μ^2 / x1 : x1
 end
 
-function rand!(rng::AbstractRNG, d::InverseGaussian, A::AbstractArray{<:Real})
-    # Based off of the non-vectorized code
-    μ, λ = params(d)
-    Z = V = W = X1 = A # prevents extra heap allocs
-    randn!(rng, A)
-    V .*= Z
-    W .*= μ
-    X1 .= μ .+ μ / (2λ) .* (W .- sqrt.(W .* (4λ .+ W)))
-    A .= ifelse.( # TODO Can we avoid heap-allocing a whole float vec?
-        rand(rng, size(A)...) .>= μ ./ (μ .+ X1),
-        (μ * μ) ./ X1,
-        X1
-    )
-end
-
-
 #### Fit model
 """
 Sufficient statistics for `InverseGaussian`, containing the weighted
diff --git a/src/univariate/continuous/laplace.jl b/src/univariate/continuous/laplace.jl
index 5bb6417ad7..e58cf6bfa0 100644
--- a/src/univariate/continuous/laplace.jl
+++ b/src/univariate/continuous/laplace.jl
@@ -126,15 +126,6 @@ Base.:*(c::Real, d::Laplace) = Laplace(c * d.μ, abs(c) * d.θ)
 rand(rng::AbstractRNG, d::Laplace) =
     d.μ + d.θ*randexp(rng)*ifelse(rand(rng, Bool), 1, -1)
 
-function rand!(rng::AbstractRNG, d::Laplace, A::AbstractArray{<:Real})
-    randexp!(rng, A)
-    A .= muladd.(
-        A,
-        ifelse.(rand(rng, Bool, size(A)), d.θ, -d.θ),
-        d.μ
-    )
-end
-
 #### Fitting
 
 function fit_mle(::Type{<:Laplace}, x::AbstractArray{<:Real})
diff --git a/src/univariate/continuous/lognormal.jl b/src/univariate/continuous/lognormal.jl
index 6969b4c9d8..afea89fb12 100644
--- a/src/univariate/continuous/lognormal.jl
+++ b/src/univariate/continuous/lognormal.jl
@@ -156,7 +156,7 @@ end
 
 #### Sampling
 
-xval(d::LogNormal, z::Real) = exp(muladd(d.σ, z, d.μ)  
+xval(d::LogNormal, z::Real) = exp(muladd(d.σ, z, d.μ))
 
 rand(rng::AbstractRNG, d::LogNormal) = xval(d,  randn(rng))
 function rand!(rng::AbstractRNG, d::LogNormal, A::AbstractArray{<:Real})
diff --git a/src/univariate/continuous/pgeneralizedgaussian.jl b/src/univariate/continuous/pgeneralizedgaussian.jl
index 4c9dcfa5ff..3c9c18edef 100644
--- a/src/univariate/continuous/pgeneralizedgaussian.jl
+++ b/src/univariate/continuous/pgeneralizedgaussian.jl
@@ -146,15 +146,4 @@ function rand(rng::AbstractRNG, d::PGeneralizedGaussian)
     else
         return d.μ + z
     end
-end
-
-function rand!(rng::AbstractRNG, d::PGeneralizedGaussian, A::AbstractArray{<:Real})
-    inv_p = inv(d.p)
-    g = Gamma(inv_p, 1)
-    A .= rand!(rng, g, A).^inv_p
-    A .= muladd.(
-        A,
-        ifelse.(rand(rng, Bool, size(A)), d.α, -d.α),
-        d.μ
-    )
 end
\ No newline at end of file

From fe2b25e5b62389d09e58abf5df29ab8009947d09 Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Sat, 17 Aug 2024 20:15:03 -0400
Subject: [PATCH 07/16] Remove redundant new tests

---
 test/univariate/continuous/inversegaussian.jl      | 14 --------------
 test/univariate/continuous/pgeneralizedgaussian.jl |  3 ---
 2 files changed, 17 deletions(-)

diff --git a/test/univariate/continuous/inversegaussian.jl b/test/univariate/continuous/inversegaussian.jl
index 972c71c2ee..6990cae710 100644
--- a/test/univariate/continuous/inversegaussian.jl
+++ b/test/univariate/continuous/inversegaussian.jl
@@ -1,17 +1,3 @@
-# Sampling Tests
-@testset "InverseGaussian sampling tests" begin
-    for d in [
-        InverseGaussian()
-        InverseGaussian(0.8)
-        InverseGaussian(2.0)
-        InverseGaussian(1.0, 1.0)
-        InverseGaussian(2.0, 1.5)
-        InverseGaussian(2.0, 7.0)
-    ]
-        test_distr(d, 10^6, test_scalar_rand = true)
-    end
-end
-
 @testset "InverseGaussian cdf outside of [0, 1] (#1873)" begin
     for d in [
         InverseGaussian(1.65, 590),
diff --git a/test/univariate/continuous/pgeneralizedgaussian.jl b/test/univariate/continuous/pgeneralizedgaussian.jl
index 52ba0ee6b0..69751155a9 100644
--- a/test/univariate/continuous/pgeneralizedgaussian.jl
+++ b/test/univariate/continuous/pgeneralizedgaussian.jl
@@ -102,8 +102,5 @@ using Test
         @test cdf(d, Inf) == 1
         @test logcdf(d, Inf) == 0
         @test quantile(d, 1 // 2) ≈ μ
-
-        # Additional tests, including sampling
-        test_distr(d, 10^6, test_scalar_rand = true)
     end
 end

From 3cd853e698c7b8703ce208a8a6320ec77bae812d Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Sat, 17 Aug 2024 20:36:32 -0400
Subject: [PATCH 08/16] Clean up more

---
 src/univariate/continuous/inversegaussian.jl       | 1 +
 src/univariate/continuous/laplace.jl               | 1 +
 test/univariate/continuous/pgeneralizedgaussian.jl | 2 +-
 3 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/univariate/continuous/inversegaussian.jl b/src/univariate/continuous/inversegaussian.jl
index ec507bc610..b585ec3b3d 100644
--- a/src/univariate/continuous/inversegaussian.jl
+++ b/src/univariate/continuous/inversegaussian.jl
@@ -170,6 +170,7 @@ function rand(rng::AbstractRNG, d::InverseGaussian)
 end
 
 #### Fit model
+
 """
 Sufficient statistics for `InverseGaussian`, containing the weighted
 sum of observations, the weighted sum of inverse points and sum of weights.
diff --git a/src/univariate/continuous/laplace.jl b/src/univariate/continuous/laplace.jl
index e58cf6bfa0..2a7bf04a47 100644
--- a/src/univariate/continuous/laplace.jl
+++ b/src/univariate/continuous/laplace.jl
@@ -126,6 +126,7 @@ Base.:*(c::Real, d::Laplace) = Laplace(c * d.μ, abs(c) * d.θ)
 rand(rng::AbstractRNG, d::Laplace) =
     d.μ + d.θ*randexp(rng)*ifelse(rand(rng, Bool), 1, -1)
 
+
 #### Fitting
 
 function fit_mle(::Type{<:Laplace}, x::AbstractArray{<:Real})
diff --git a/test/univariate/continuous/pgeneralizedgaussian.jl b/test/univariate/continuous/pgeneralizedgaussian.jl
index 69751155a9..9857e4fb9a 100644
--- a/test/univariate/continuous/pgeneralizedgaussian.jl
+++ b/test/univariate/continuous/pgeneralizedgaussian.jl
@@ -74,7 +74,7 @@ using Test
             end
 
             # Additional tests, including sampling
-            test_distr(d, 10^6, test_scalar_rand = true)
+            test_distr(d, 10^6)
         end
     end
 

From d53507fc5cb6798a6b66cd098248a6cec178be89 Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Sat, 17 Aug 2024 20:38:12 -0400
Subject: [PATCH 09/16] Partially undo previous undo to changes to tests

---
 test/univariate/continuous/pgeneralizedgaussian.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/test/univariate/continuous/pgeneralizedgaussian.jl b/test/univariate/continuous/pgeneralizedgaussian.jl
index 9857e4fb9a..49de167cc5 100644
--- a/test/univariate/continuous/pgeneralizedgaussian.jl
+++ b/test/univariate/continuous/pgeneralizedgaussian.jl
@@ -102,5 +102,8 @@ using Test
         @test cdf(d, Inf) == 1
         @test logcdf(d, Inf) == 0
         @test quantile(d, 1 // 2) ≈ μ
+
+        # Additional tests, including sampling
+        test_distr(d, 10^6)
     end
 end

From efc0a60fa70fef271326db052bd324219735a347 Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Mon, 26 Aug 2024 21:27:54 -0400
Subject: [PATCH 10/16] Use xval for NormalCanon rand

---
 src/univariate/continuous/normalcanon.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/univariate/continuous/normalcanon.jl b/src/univariate/continuous/normalcanon.jl
index 36f1215be2..89d45a4fa6 100644
--- a/src/univariate/continuous/normalcanon.jl
+++ b/src/univariate/continuous/normalcanon.jl
@@ -87,11 +87,11 @@ invlogccdf(d::NormalCanon, lp::Real) = xval(d, norminvlogccdf(lp))
 
 #### Sampling
 
-rand(rng::AbstractRNG, cf::NormalCanon) = cf.μ + randn(rng) / sqrt(cf.λ)
+rand(rng::AbstractRNG, cf::NormalCanon) = xval(cf, randn(rng))
 
 function rand!(rng::AbstractRNG, cf::NormalCanon, A::AbstractArray{<:Real})
     randn!(rng, A)
-    A .= cf.μ .+ A ./ sqrt(cf.λ)
+    map!(Base.Fix1(xval, cf), A, A)
     return A
 end
 

From d13c8fb7604dcefa83693e38594718a2aa4f4285 Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Tue, 3 Sep 2024 23:28:33 -0400
Subject: [PATCH 11/16] Apply suggestions from code review

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
---
 test/testutils.jl                         | 19 +++++++------------
 test/univariate/continuous/exponential.jl |  2 +-
 test/univariate/continuous/logitnormal.jl |  2 +-
 test/univariate/continuous/lognormal.jl   |  2 +-
 test/univariate/continuous/normalcanon.jl |  2 +-
 test/univariate/continuous/pareto.jl      |  2 +-
 6 files changed, 12 insertions(+), 17 deletions(-)

diff --git a/test/testutils.jl b/test/testutils.jl
index 358eb6c952..3aef74cd2b 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -30,7 +30,6 @@ end
 function test_distr(distr::DiscreteUnivariateDistribution, n::Int;
                     testquan::Bool=true, rng::AbstractRNG=MersenneTwister(),
                     test_scalar_rand::Bool=false)
-    println("    testing $(distr)")
     test_range(distr)
     vs = get_evalsamples(distr, 0.00001)
 
@@ -73,9 +72,7 @@ end
 # testing the implementation of a continuous univariate distribution
 #
 function test_distr(distr::ContinuousUnivariateDistribution, n::Int;
-                    testquan::Bool=true, rng::AbstractRNG=MersenneTwister(123),
-                    test_scalar_rand::Bool=false)
-    println("    testing $(distr)")
+                    testquan::Bool=true, rng::AbstractRNG=MersenneTwister(123))
     test_range(distr)
     vs = get_evalsamples(distr, 0.01, 2000)
 
@@ -116,8 +113,7 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
                       n::Int;                                   # number of samples to generate
                       q::Float64=1.0e-7,                        # confidence interval, 1 - q as confidence
                       verbose::Bool=false,                      # show intermediate info (for debugging)
-                      rng::Union{AbstractRNG, Missing}=missing,  # add an rng?
-                      call_scalar::Bool=false)                  # directly scall the scalar rand(d) instead
+                      rng::Union{AbstractRNG, Missing}=missing) # add an rng?
     # The basic idea
     # ------------------
     #   Generate n samples, and count the occurrences of each value within a reasonable range.
@@ -211,8 +207,8 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
 end
 
 test_samples(distr::DiscreteUnivariateDistribution, n::Int;
-             q::Float64=1.0e-6, verbose::Bool=false, rng=missing, kwargs...) =
-    test_samples(distr, distr, n; q=q, verbose=verbose, rng=rng, kwargs...)
+             q::Float64=1.0e-6, verbose::Bool=false, rng=missing) =
+    test_samples(distr, distr, n; q=q, verbose=verbose, rng=rng)
 
 # for continuous samplers
 #
@@ -222,8 +218,7 @@ function test_samples(s::Sampleable{Univariate, Continuous},    # the sampleable
                       nbins::Int=50,                            # divide the main interval into nbins
                       q::Float64=1.0e-6,                        # confidence interval, 1 - q as confidence
                       verbose::Bool=false,                      # show intermediate info (for debugging)
-                      rng::Union{AbstractRNG, Missing}=missing, # add an rng?
-                      call_scalar::Bool=false)                  # directly scall the scalar rand(d) instead
+                      rng::Union{AbstractRNG, Missing}=missing) # add an rng?
 
     # The basic idea
     # ------------------
@@ -332,8 +327,8 @@ function test_samples(s::Sampleable{Univariate, Continuous},    # the sampleable
     return samples
 end
 
-test_samples(distr::ContinuousUnivariateDistribution, n::Int; nbins::Int=50, q::Float64=1.0e-6, verbose::Bool=false, rng=missing, kwargs...) =
-    test_samples(distr, distr, n; nbins=nbins, q=q, verbose=verbose, rng=rng, kwargs...)
+test_samples(distr::ContinuousUnivariateDistribution, n::Int; nbins::Int=50, q::Float64=1.0e-6, verbose::Bool=false, rng=missing) =
+    test_samples(distr, distr, n; nbins=nbins, q=q, verbose=verbose, rng=rng)
 
 
 #### Testing range & support methods
diff --git a/test/univariate/continuous/exponential.jl b/test/univariate/continuous/exponential.jl
index e8284389a1..e710dadaba 100644
--- a/test/univariate/continuous/exponential.jl
+++ b/test/univariate/continuous/exponential.jl
@@ -19,6 +19,6 @@ test_cgf(Exponential(10  ), (0.08, -1, -100f0, -1e6))
         Exponential(0.91),
         Exponential(10)
     ]
-        test_distr(d, 10^6, test_scalar_rand = true)
+        test_distr(d, 10^6)
     end
 end
diff --git a/test/univariate/continuous/logitnormal.jl b/test/univariate/continuous/logitnormal.jl
index b8c4a578d8..9959608d33 100644
--- a/test/univariate/continuous/logitnormal.jl
+++ b/test/univariate/continuous/logitnormal.jl
@@ -70,6 +70,6 @@ end
         LogitNormal(-2, 3),
         LogitNormal(0, 0.2)
     ]
-        test_distr(d, 10^6, test_scalar_rand = true)
+        test_distr(d, 10^6)
     end
 end
\ No newline at end of file
diff --git a/test/univariate/continuous/lognormal.jl b/test/univariate/continuous/lognormal.jl
index 26c16f7730..3564d6fb87 100644
--- a/test/univariate/continuous/lognormal.jl
+++ b/test/univariate/continuous/lognormal.jl
@@ -325,6 +325,6 @@ end
             LogNormal(3.0, 1.0)
             LogNormal(3.0, 2.0)
         ]
-        test_distr(d, 10^6, test_scalar_rand = true)
+        test_distr(d, 10^6)
     end
 end
\ No newline at end of file
diff --git a/test/univariate/continuous/normalcanon.jl b/test/univariate/continuous/normalcanon.jl
index 443d6ec494..4fc8de4daa 100644
--- a/test/univariate/continuous/normalcanon.jl
+++ b/test/univariate/continuous/normalcanon.jl
@@ -5,6 +5,6 @@
         NormalCanon(-1.0, 2.5)
         NormalCanon(2.0, 0.8)
     ]
-        test_distr(d, 10^6, test_scalar_rand = true)
+        test_distr(d, 10^6)
     end
 end
\ No newline at end of file
diff --git a/test/univariate/continuous/pareto.jl b/test/univariate/continuous/pareto.jl
index 9e7e95775e..87ef527c7a 100644
--- a/test/univariate/continuous/pareto.jl
+++ b/test/univariate/continuous/pareto.jl
@@ -5,6 +5,6 @@
             Pareto(2.0, 1.5)
             Pareto(3.0, 2.0)
         ]
-        test_distr(d, 10^6, test_scalar_rand = true)
+        test_distr(d, 10^6)
     end
 end
\ No newline at end of file

From 1f7e28a9736a17167f56c4938402d3f38caa989d Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Tue, 3 Sep 2024 23:58:04 -0400
Subject: [PATCH 12/16] Apply other recommendations to testutils

---
 test/testutils.jl | 113 ++++++++++++++++++++++------------------------
 1 file changed, 55 insertions(+), 58 deletions(-)

diff --git a/test/testutils.jl b/test/testutils.jl
index 3aef74cd2b..b922c312cb 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -28,8 +28,7 @@ end
 # testing the implementation of a discrete univariate distribution
 #
 function test_distr(distr::DiscreteUnivariateDistribution, n::Int;
-                    testquan::Bool=true, rng::AbstractRNG=MersenneTwister(),
-                    test_scalar_rand::Bool=false)
+                    testquan::Bool=true, rng::AbstractRNG=MersenneTwister(123))
     test_range(distr)
     vs = get_evalsamples(distr, 0.00001)
 
@@ -41,10 +40,6 @@ function test_distr(distr::DiscreteUnivariateDistribution, n::Int;
     test_stats(distr, vs)
     test_samples(distr, n)
     test_samples(distr, n, rng=rng)
-    if test_scalar_rand
-        xs = test_samples(distr, n; call_scalar = true)
-        xs = test_samples(distr, n, rng=rng, call_scalar = true)
-    end
 
     test_params(distr)
 end
@@ -87,12 +82,6 @@ function test_distr(distr::ContinuousUnivariateDistribution, n::Int;
     allow_test_stats(distr) && test_stats(distr, xs)
     xs = test_samples(distr, n, rng=rng)
     allow_test_stats(distr) && test_stats(distr, xs)
-    if test_scalar_rand
-        xs = test_samples(distr, n; call_scalar = true)
-        allow_test_stats(distr) && test_stats(distr, xs)
-        xs = test_samples(distr, n, rng=rng, call_scalar = true)
-        allow_test_stats(distr) && test_stats(distr, xs)
-    end
 
     test_params(distr)
 end
@@ -157,43 +146,45 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
     # generate samples using RNG passed or default RNG
     # we also check reproducibility
     if rng === missing
-        Random.seed!(1234)    
-        samples = if !call_scalar
-            rand(s, n)
-        else
-            map((_) -> rand(s), 1:n)
-        end
         Random.seed!(1234)
-        samples2 = if !call_scalar
-            rand(s, n)
-        else
-            map((_) -> rand(s), 1:n)
-        end
+        samples = rand(s, n)
+        Random.seed!(1234)
+        samples2 = rand(s, n)
+        Random.seed!(1234)
+        samples3 = map!((_) -> rand(s), 1:n)
+        Random.seed!(1234)
+        samples4 = map((_) -> rand(s), 1:n)
     else
         rng2 = deepcopy(rng)
-        samples = if !call_scalar
-            rand(rng, s, n)
-        else
-            map((_) -> rand(rng, s), 1:n)
-        end
-        samples2 = if !call_scalar
-            rand(rng2, s, n)
-        else
-            map((_) -> rand(rng2, s), 1:n)
-        end
+        rng3 = deepcopy(rng)
+        rng4 = deepcopy(rng)
+        samples = rand(rng, s, n)
+        samples2 = rand(rng2, s, n)
+        samples3 = map((_) -> rand(rng3, s), 1:n)
+        samples4 = map((_) -> rand(rng4, s), 1:n)
     end
     @test length(samples) == n
     @test samples2 == samples
+    @test samples3 == samples4
 
     # scan samples and get counts
     cnts = zeros(Int, m)
+    cnts_sc = zeros(Int, m)
     for i = 1:n
         @inbounds si = samples[i]
         if rmin <= si <= rmax
             cnts[si - rmin + 1] += 1
         else
             vmin <= si <= vmax ||
-                error("Sample value out of valid range.")
+                error("Sample value out of valid range. (Vector Method)")
+        end
+
+        @inbounds si_sc = samples3[i]
+        if rmin <= si_sc <= rmax
+            cnts_sc[si_sc - rmin + 1] += 1
+        else
+            vmin <= si_sc <= vmax ||
+                error("Sample value out of valid range. (Scalar Method)")
         end
     end
 
@@ -201,7 +192,11 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
     for i = 1:m
         verbose && println("v = $(rmin+i-1) ==> ($(clb[i]), $(cub[i])): $(cnts[i])")
         clb[i] <= cnts[i] <= cub[i] ||
-            error("The counts are out of the confidence interval.")
+            error("The counts are out of the confidence interval. (Vector Method)")
+
+        verbose && println("v = $(rmin+i-1) ==> ($(clb[i]), $(cub[i])): $(cnts_sc[i])")
+        clb[i] <= cnts_sc[i] <= cub[i] ||
+            error("The counts are out of the confidence interval. (Scalar Method)")
     end
     return samples
 end
@@ -273,33 +268,26 @@ function test_samples(s::Sampleable{Univariate, Continuous},    # the sampleable
     # generate samples using RNG passed or default RNG
     # we also check reproducibility
     if rng === missing
-        Random.seed!(1234)    
-        samples = if !call_scalar
-            rand(s, n)
-        else
-            map((_) -> rand(s), 1:n)
-        end
         Random.seed!(1234)
-        samples2 = if !call_scalar
-            rand(s, n)
-        else
-            map((_) -> rand(s), 1:n)
-        end
+        samples = rand(s, n)
+        Random.seed!(1234)
+        samples2 = rand(s, n)
+        Random.seed!(1234)
+        samples3 = map!((_) -> rand(s), 1:n)
+        Random.seed!(1234)
+        samples4 = map((_) -> rand(s), 1:n)
     else
         rng2 = deepcopy(rng)
-        samples = if !call_scalar
-            rand(rng, s, n)
-        else
-            map((_) -> rand(rng, s), 1:n)
-        end
-        samples2 = if !call_scalar
-            rand(rng2, s, n)
-        else
-            map((_) -> rand(rng2, s), 1:n)
-        end
+        rng3 = deepcopy(rng)
+        rng4 = deepcopy(rng)
+        samples = rand(rng, s, n)
+        samples2 = rand(rng2, s, n)
+        samples3 = map((_) -> rand(rng3, s), 1:n)
+        samples4 = map((_) -> rand(rng4, s), 1:n)
     end
     @test length(samples) == n
     @test samples2 == samples
+    @test samples3 == samples4
 
     if isa(distr, StudentizedRange)
         samples[isnan.(samples)] .= 0.0 # Underlying implementation in Rmath can't handle very low values.
@@ -309,20 +297,29 @@ function test_samples(s::Sampleable{Univariate, Continuous},    # the sampleable
     for i = 1:n
         @inbounds si = samples[i]
         vmin <= si <= vmax ||
-            error("Sample value out of valid range.")
+            error("Sample value out of valid range. (Vector Method)")
+        @inbounds si_sc = samples3[i]
+        vmin <= si_sc <= vmax ||
+            error("Sample value out of valid range. (Scalar Method)")
     end
 
     # get counts
     cnts = fit(Histogram, samples, edges; closed=:right).weights
     @assert length(cnts) == nbins
 
+    cnts_sc = fit(Histogram, samples3, edges; closed=:right).weights
+    @assert length(cnts_sc) == nbins
+
     # check the counts
     for i = 1:nbins
         if verbose
             @printf("[%.4f, %.4f) ==> (%d, %d): %d\n", edges[i], edges[i+1], clb[i], cub[i], cnts[i])
+            @printf("[%.4f, %.4f) ==> (%d, %d): %d\n", edges[i], edges[i+1], clb[i], cub[i], cnts_sc[i])
         end
         clb[i] <= cnts[i] <= cub[i] ||
-            error("The counts are out of the confidence interval.")
+            error("The counts are out of the confidence interval. (Vector Method)")
+        clb[i] <= cnts_sc[i] <= cub[i] ||
+            error("The counts are out of the confidence interval. (Scalar Method)")
     end
     return samples
 end

From 0746991c252281f07af112b79a9ba5bdab08c3f1 Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Wed, 4 Sep 2024 00:22:53 -0400
Subject: [PATCH 13/16] Fix erroneous !

---
 test/testutils.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/testutils.jl b/test/testutils.jl
index b922c312cb..7e72679966 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -151,7 +151,7 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
         Random.seed!(1234)
         samples2 = rand(s, n)
         Random.seed!(1234)
-        samples3 = map!((_) -> rand(s), 1:n)
+        samples3 = map((_) -> rand(s), 1:n)
         Random.seed!(1234)
         samples4 = map((_) -> rand(s), 1:n)
     else
@@ -273,7 +273,7 @@ function test_samples(s::Sampleable{Univariate, Continuous},    # the sampleable
         Random.seed!(1234)
         samples2 = rand(s, n)
         Random.seed!(1234)
-        samples3 = map!((_) -> rand(s), 1:n)
+        samples3 = map((_) -> rand(s), 1:n)
         Random.seed!(1234)
         samples4 = map((_) -> rand(s), 1:n)
     else

From 38faf5ca78778ffb22b4587f662be0286d9fa23c Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Wed, 25 Sep 2024 10:35:13 +0200
Subject: [PATCH 14/16] Address reviewer comments

---
 src/univariate/continuous/lognormal.jl        |  2 +-
 .../continuous/pgeneralizedgaussian.jl        |  2 +-
 test/fit.jl                                   |  2 +-
 test/testutils.jl                             | 41 +++++++++----------
 test/univariate/continuous/exponential.jl     |  2 +-
 test/univariate/continuous/logitnormal.jl     |  2 +-
 test/univariate/continuous/lognormal.jl       |  2 +-
 test/univariate/continuous/normalcanon.jl     |  2 +-
 test/univariate/continuous/pareto.jl          |  2 +-
 9 files changed, 28 insertions(+), 29 deletions(-)

diff --git a/src/univariate/continuous/lognormal.jl b/src/univariate/continuous/lognormal.jl
index afea89fb12..d9ada052f6 100644
--- a/src/univariate/continuous/lognormal.jl
+++ b/src/univariate/continuous/lognormal.jl
@@ -158,7 +158,7 @@ end
 
 xval(d::LogNormal, z::Real) = exp(muladd(d.σ, z, d.μ))
 
-rand(rng::AbstractRNG, d::LogNormal) = xval(d,  randn(rng))
+rand(rng::AbstractRNG, d::LogNormal) = xval(d, randn(rng))
 function rand!(rng::AbstractRNG, d::LogNormal, A::AbstractArray{<:Real})
     randn!(rng, A)
     map!(Base.Fix1(xval, d), A, A)
diff --git a/src/univariate/continuous/pgeneralizedgaussian.jl b/src/univariate/continuous/pgeneralizedgaussian.jl
index 3c9c18edef..24eb9a9698 100644
--- a/src/univariate/continuous/pgeneralizedgaussian.jl
+++ b/src/univariate/continuous/pgeneralizedgaussian.jl
@@ -146,4 +146,4 @@ function rand(rng::AbstractRNG, d::PGeneralizedGaussian)
     else
         return d.μ + z
     end
-end
\ No newline at end of file
+end
diff --git a/test/fit.jl b/test/fit.jl
index b3474c6fe8..143f0b6ee4 100644
--- a/test/fit.jl
+++ b/test/fit.jl
@@ -369,7 +369,7 @@ end
     for func in funcs, dist in (Laplace, Laplace{Float64})
         d = fit(dist, func[2](dist(5.0, 3.0), N + 1))
         @test isa(d, dist)
-        @test isapprox(location(d), 5.0, atol=0.021)
+        @test isapprox(location(d), 5.0, atol=0.03)
         @test isapprox(scale(d)   , 3.0, atol=0.03)
     end
 end
diff --git a/test/testutils.jl b/test/testutils.jl
index 7e72679966..140eaab729 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -28,7 +28,8 @@ end
 # testing the implementation of a discrete univariate distribution
 #
 function test_distr(distr::DiscreteUnivariateDistribution, n::Int;
-                    testquan::Bool=true, rng::AbstractRNG=MersenneTwister(123))
+                    testquan::Bool=true, rng::AbstractRNG = Random.default_rng())
+
     test_range(distr)
     vs = get_evalsamples(distr, 0.00001)
 
@@ -39,8 +40,7 @@ function test_distr(distr::DiscreteUnivariateDistribution, n::Int;
 
     test_stats(distr, vs)
     test_samples(distr, n)
-    test_samples(distr, n, rng=rng)
-
+    test_samples(distr, n; rng=rng)
     test_params(distr)
 end
 
@@ -82,7 +82,6 @@ function test_distr(distr::ContinuousUnivariateDistribution, n::Int;
     allow_test_stats(distr) && test_stats(distr, xs)
     xs = test_samples(distr, n, rng=rng)
     allow_test_stats(distr) && test_stats(distr, xs)
-
     test_params(distr)
 end
 
@@ -103,6 +102,7 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
                       q::Float64=1.0e-7,                        # confidence interval, 1 - q as confidence
                       verbose::Bool=false,                      # show intermediate info (for debugging)
                       rng::Union{AbstractRNG, Missing}=missing) # add an rng?
+
     # The basic idea
     # ------------------
     #   Generate n samples, and count the occurrences of each value within a reasonable range.
@@ -151,17 +151,17 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
         Random.seed!(1234)
         samples2 = rand(s, n)
         Random.seed!(1234)
-        samples3 = map((_) -> rand(s), 1:n)
+        samples3 = [rand(s) for _ in 1:n]
         Random.seed!(1234)
-        samples4 = map((_) -> rand(s), 1:n)
+        samples4 = [rand(s) for _ in 1:n]
     else
         rng2 = deepcopy(rng)
         rng3 = deepcopy(rng)
         rng4 = deepcopy(rng)
         samples = rand(rng, s, n)
         samples2 = rand(rng2, s, n)
-        samples3 = map((_) -> rand(rng3, s), 1:n)
-        samples4 = map((_) -> rand(rng4, s), 1:n)
+        samples3 = [rand(rng3, s) for _ in 1:n]
+        samples4 = [rand(rng4, s) for _ in 1:n]
     end
     @test length(samples) == n
     @test samples2 == samples
@@ -176,7 +176,7 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
             cnts[si - rmin + 1] += 1
         else
             vmin <= si <= vmax ||
-                error("Sample value out of valid range. (Vector Method)")
+                throw(DomainError(si, "sample generated by `rand(s, n)` is out of valid range [$vmin, $vmax]."))
         end
 
         @inbounds si_sc = samples3[i]
@@ -184,7 +184,7 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
             cnts_sc[si_sc - rmin + 1] += 1
         else
             vmin <= si_sc <= vmax ||
-                error("Sample value out of valid range. (Scalar Method)")
+                throw(DomainError(si, "sample generated by `[rand(s) for _ in 1:n]` is out of valid range [$vmin, $vmax]."))
         end
     end
 
@@ -192,11 +192,11 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
     for i = 1:m
         verbose && println("v = $(rmin+i-1) ==> ($(clb[i]), $(cub[i])): $(cnts[i])")
         clb[i] <= cnts[i] <= cub[i] ||
-            error("The counts are out of the confidence interval. (Vector Method)")
+            error("The counts of samples generated by `rand(s, n)` are out of the confidence interval.")
 
         verbose && println("v = $(rmin+i-1) ==> ($(clb[i]), $(cub[i])): $(cnts_sc[i])")
         clb[i] <= cnts_sc[i] <= cub[i] ||
-            error("The counts are out of the confidence interval. (Scalar Method)")
+            error("The counts of samples generated by `[rand(s) for _ in 1:n]` are out of the confidence interval.")
     end
     return samples
 end
@@ -273,17 +273,17 @@ function test_samples(s::Sampleable{Univariate, Continuous},    # the sampleable
         Random.seed!(1234)
         samples2 = rand(s, n)
         Random.seed!(1234)
-        samples3 = map((_) -> rand(s), 1:n)
+        samples3 = [rand(s) for _ in 1:n]
         Random.seed!(1234)
-        samples4 = map((_) -> rand(s), 1:n)
+        samples4 = [rand(s) for _ in 1:n]
     else
         rng2 = deepcopy(rng)
         rng3 = deepcopy(rng)
         rng4 = deepcopy(rng)
         samples = rand(rng, s, n)
         samples2 = rand(rng2, s, n)
-        samples3 = map((_) -> rand(rng3, s), 1:n)
-        samples4 = map((_) -> rand(rng4, s), 1:n)
+        samples3 = [rand(rng3, s) for _ in 1:n]
+        samples4 = [rand(rng4, s) for _ in 1:n]
     end
     @test length(samples) == n
     @test samples2 == samples
@@ -297,10 +297,10 @@ function test_samples(s::Sampleable{Univariate, Continuous},    # the sampleable
     for i = 1:n
         @inbounds si = samples[i]
         vmin <= si <= vmax ||
-            error("Sample value out of valid range. (Vector Method)")
+            throw(DomainError(si, "sample generated by `rand(s, n)` is out of valid range [$vmin, $vmax]."))
         @inbounds si_sc = samples3[i]
         vmin <= si_sc <= vmax ||
-            error("Sample value out of valid range. (Scalar Method)")
+            throw(DomainError(si, "sample generated by `[rand(s) for _ in 1:n]` is out of valid range [$vmin, $vmax]."))
     end
 
     # get counts
@@ -317,9 +317,9 @@ function test_samples(s::Sampleable{Univariate, Continuous},    # the sampleable
             @printf("[%.4f, %.4f) ==> (%d, %d): %d\n", edges[i], edges[i+1], clb[i], cub[i], cnts_sc[i])
         end
         clb[i] <= cnts[i] <= cub[i] ||
-            error("The counts are out of the confidence interval. (Vector Method)")
+            error("The counts of samples generated by `rand(s, n)` are out of the confidence interval.")
         clb[i] <= cnts_sc[i] <= cub[i] ||
-            error("The counts are out of the confidence interval. (Scalar Method)")
+            error("The counts of samples generated by `[rand(s) for _ in 1:n]` are out of the confidence interval.")
     end
     return samples
 end
@@ -623,7 +623,6 @@ end
 allow_test_stats(d::UnivariateDistribution) = true
 allow_test_stats(d::NoncentralBeta) = false
 allow_test_stats(::StudentizedRange) = false
-allow_test_stats(::LogitNormal) = false
 
 function test_stats(d::ContinuousUnivariateDistribution, xs::AbstractVector{Float64})
     # using Monte Carlo methods
diff --git a/test/univariate/continuous/exponential.jl b/test/univariate/continuous/exponential.jl
index e710dadaba..9c5f29aa93 100644
--- a/test/univariate/continuous/exponential.jl
+++ b/test/univariate/continuous/exponential.jl
@@ -10,7 +10,7 @@ end
 
 test_cgf(Exponential(1), (0.9, -1, -100f0, -1e6))
 test_cgf(Exponential(0.91), (0.9, -1, -100f0, -1e6))
-test_cgf(Exponential(10  ), (0.08, -1, -100f0, -1e6))
+test_cgf(Exponential(10), (0.08, -1, -100f0, -1e6))
 
 # Sampling Tests
 @testset "Exponential sampling tests" begin
diff --git a/test/univariate/continuous/logitnormal.jl b/test/univariate/continuous/logitnormal.jl
index 9959608d33..b8d72290ef 100644
--- a/test/univariate/continuous/logitnormal.jl
+++ b/test/univariate/continuous/logitnormal.jl
@@ -72,4 +72,4 @@ end
     ]
         test_distr(d, 10^6)
     end
-end
\ No newline at end of file
+end
diff --git a/test/univariate/continuous/lognormal.jl b/test/univariate/continuous/lognormal.jl
index 3564d6fb87..1cc44b433d 100644
--- a/test/univariate/continuous/lognormal.jl
+++ b/test/univariate/continuous/lognormal.jl
@@ -327,4 +327,4 @@ end
         ]
         test_distr(d, 10^6)
     end
-end
\ No newline at end of file
+end
diff --git a/test/univariate/continuous/normalcanon.jl b/test/univariate/continuous/normalcanon.jl
index 4fc8de4daa..9ae52e15f0 100644
--- a/test/univariate/continuous/normalcanon.jl
+++ b/test/univariate/continuous/normalcanon.jl
@@ -7,4 +7,4 @@
     ]
         test_distr(d, 10^6)
     end
-end
\ No newline at end of file
+end
diff --git a/test/univariate/continuous/pareto.jl b/test/univariate/continuous/pareto.jl
index 87ef527c7a..195e34a96e 100644
--- a/test/univariate/continuous/pareto.jl
+++ b/test/univariate/continuous/pareto.jl
@@ -7,4 +7,4 @@
         ]
         test_distr(d, 10^6)
     end
-end
\ No newline at end of file
+end

From 3292a5acd0e482f6d866334a0b473bd4b2b2fa70 Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Wed, 25 Sep 2024 11:27:04 +0200
Subject: [PATCH 15/16] `mean` not defined for `LogitNormal`

---
 test/testutils.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/test/testutils.jl b/test/testutils.jl
index 140eaab729..85d2c6863b 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -623,6 +623,7 @@ end
 allow_test_stats(d::UnivariateDistribution) = true
 allow_test_stats(d::NoncentralBeta) = false
 allow_test_stats(::StudentizedRange) = false
+allow_test_stats(::LogitNormal) = false # `mean` is not defined since it has no analytical solution
 
 function test_stats(d::ContinuousUnivariateDistribution, xs::AbstractVector{Float64})
     # using Monte Carlo methods

From a704c6e71ed483f9cdce5b026f63433f6af49efa Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Wed, 25 Sep 2024 11:59:24 +0200
Subject: [PATCH 16/16] Copy RNG with `copy`, not `deepcopy`

---
 test/testutils.jl | 16 ++++++++++------
 1 file changed, 10 insertions(+), 6 deletions(-)

diff --git a/test/testutils.jl b/test/testutils.jl
index 85d2c6863b..8acd378535 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -155,9 +155,11 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
         Random.seed!(1234)
         samples4 = [rand(s) for _ in 1:n]
     else
-        rng2 = deepcopy(rng)
-        rng3 = deepcopy(rng)
-        rng4 = deepcopy(rng)
+        # RNGs have to be copied with `copy`, not `deepcopy`
+        # Ref https://github.com/JuliaLang/julia/issues/42899
+        rng2 = copy(rng)
+        rng3 = copy(rng)
+        rng4 = copy(rng)
         samples = rand(rng, s, n)
         samples2 = rand(rng2, s, n)
         samples3 = [rand(rng3, s) for _ in 1:n]
@@ -277,9 +279,11 @@ function test_samples(s::Sampleable{Univariate, Continuous},    # the sampleable
         Random.seed!(1234)
         samples4 = [rand(s) for _ in 1:n]
     else
-        rng2 = deepcopy(rng)
-        rng3 = deepcopy(rng)
-        rng4 = deepcopy(rng)
+        # RNGs have to be copied with `copy`, not `deepcopy`
+        # Ref https://github.com/JuliaLang/julia/issues/42899
+        rng2 = copy(rng)
+        rng3 = copy(rng)
+        rng4 = copy(rng)
         samples = rand(rng, s, n)
         samples2 = rand(rng2, s, n)
         samples3 = [rand(rng3, s) for _ in 1:n]