From 08c56eace1ed12bc74941e7b6c8f66bd7d3f7a3c Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Wed, 25 Sep 2024 07:46:51 -0400
Subject: [PATCH 1/9] Specialized vector rand! for many distributions (#1879)

* Test scalar rand separately from vector rand

* Add specialized rand! for many distributions

* Restore location of old NormalInverseGaussian tests

* Remove duplication of inversegaussian in runtests.jl

* Apply many suggestions from code review

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>

* Apply other suggestions

* Remove redundant new tests

* Clean up more

* Partially undo previous undo to changes to tests

* Use xval for NormalCanon rand

* Apply suggestions from code review

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>

* Apply other recommendations to testutils

* Fix erroneous !

* Address reviewer comments

* `mean` not defined for `LogitNormal`

* Copy RNG with `copy`, not `deepcopy`

---------

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
---
 src/univariate/continuous/exponential.jl      |  5 ++
 src/univariate/continuous/logitnormal.jl      |  9 ++-
 src/univariate/continuous/lognormal.jl        |  9 ++-
 src/univariate/continuous/normal.jl           |  9 ++-
 src/univariate/continuous/normalcanon.jl      |  8 ++-
 src/univariate/continuous/pareto.jl           |  9 ++-
 .../continuous/pgeneralizedgaussian.jl        |  2 +-
 test/fit.jl                                   |  2 +-
 test/multivariate/mvlognormal.jl              |  4 +-
 test/runtests.jl                              |  4 +-
 test/testutils.jl                             | 61 ++++++++++++++++---
 test/univariate/continuous/exponential.jl     | 16 ++++-
 test/univariate/continuous/logitnormal.jl     |  9 +++
 test/univariate/continuous/lognormal.jl       | 14 +++++
 test/univariate/continuous/normalcanon.jl     | 10 +++
 test/univariate/continuous/pareto.jl          | 10 +++
 16 files changed, 160 insertions(+), 21 deletions(-)
 create mode 100644 test/univariate/continuous/normalcanon.jl
 create mode 100644 test/univariate/continuous/pareto.jl

diff --git a/src/univariate/continuous/exponential.jl b/src/univariate/continuous/exponential.jl
index 573a469b07..96737c94a2 100644
--- a/src/univariate/continuous/exponential.jl
+++ b/src/univariate/continuous/exponential.jl
@@ -107,6 +107,11 @@ cf(d::Exponential, t::Real) = 1/(1 - t * im * scale(d))
 #### Sampling
 rand(rng::AbstractRNG, d::Exponential{T}) where {T} = xval(d, randexp(rng, float(T)))
 
+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 7fd5fb9113..cf0832d00c 100644
--- a/src/univariate/continuous/logitnormal.jl
+++ b/src/univariate/continuous/logitnormal.jl
@@ -157,7 +157,14 @@ end
 
 #### Sampling
 
-rand(rng::AbstractRNG, d::LogitNormal) = logistic(randn(rng) * d.σ + 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 666b14248a..d9ada052f6 100644
--- a/src/univariate/continuous/lognormal.jl
+++ b/src/univariate/continuous/lognormal.jl
@@ -156,7 +156,14 @@ end
 
 #### Sampling
 
-rand(rng::AbstractRNG, d::LogNormal) = exp(randn(rng) * d.σ + 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})
+    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 4231522674..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} = d.μ + d.σ * randn(rng, float(T))
+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 de29c4b196..89d45a4fa6 100644
--- a/src/univariate/continuous/normalcanon.jl
+++ b/src/univariate/continuous/normalcanon.jl
@@ -87,7 +87,13 @@ 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)
+    map!(Base.Fix1(xval, cf), A, A)
+    return A
+end
 
 #### Affine transformations
 
diff --git a/src/univariate/continuous/pareto.jl b/src/univariate/continuous/pareto.jl
index 10bb246c5f..249fa29eb7 100644
--- a/src/univariate/continuous/pareto.jl
+++ b/src/univariate/continuous/pareto.jl
@@ -110,7 +110,14 @@ quantile(d::Pareto, p::Real) = cquantile(d, 1 - p)
 
 #### Sampling
 
-rand(rng::AbstractRNG, d::Pareto) = d.θ * exp(randexp(rng) / 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
 
diff --git a/src/univariate/continuous/pgeneralizedgaussian.jl b/src/univariate/continuous/pgeneralizedgaussian.jl
index 9d41d9d464..24eb9a9698 100644
--- a/src/univariate/continuous/pgeneralizedgaussian.jl
+++ b/src/univariate/continuous/pgeneralizedgaussian.jl
@@ -141,7 +141,7 @@ 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
diff --git a/test/fit.jl b/test/fit.jl
index 4483dd1ec6..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.02)
+        @test isapprox(location(d), 5.0, atol=0.03)
         @test isapprox(scale(d)   , 3.0, atol=0.03)
     end
 end
diff --git a/test/multivariate/mvlognormal.jl b/test/multivariate/mvlognormal.jl
index 7ef76bef2d..af887f40d0 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 c7d73f4546..cd3afe17f9 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",
@@ -83,6 +84,7 @@ const tests = [
     "univariate/continuous/noncentralchisq",
     "univariate/continuous/weibull",
     "pdfnorm",
+    "univariate/continuous/pareto",
     "univariate/continuous/rician",
     "functionals",
     "density_interface",
@@ -143,9 +145,7 @@ const tests = [
     # "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..8acd378535 100644
--- a/test/testutils.jl
+++ b/test/testutils.jl
@@ -28,7 +28,7 @@ end
 # testing the implementation of a discrete univariate distribution
 #
 function test_distr(distr::DiscreteUnivariateDistribution, n::Int;
-                    testquan::Bool=true)
+                    testquan::Bool=true, rng::AbstractRNG = Random.default_rng())
 
     test_range(distr)
     vs = get_evalsamples(distr, 0.00001)
@@ -40,7 +40,7 @@ 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)
     test_params(distr)
 end
 
@@ -150,23 +150,43 @@ function test_samples(s::Sampleable{Univariate, Discrete},      # the sampleable
         samples = rand(s, n)
         Random.seed!(1234)
         samples2 = rand(s, n)
+        Random.seed!(1234)
+        samples3 = [rand(s) for _ in 1:n]
+        Random.seed!(1234)
+        samples4 = [rand(s) for _ in 1:n]
     else
-        rng2 = 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]
+        samples4 = [rand(rng4, s) for _ in 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.")
+                throw(DomainError(si, "sample generated by `rand(s, n)` is out of valid range [$vmin, $vmax]."))
+        end
+
+        @inbounds si_sc = samples3[i]
+        if rmin <= si_sc <= rmax
+            cnts_sc[si_sc - rmin + 1] += 1
+        else
+            vmin <= si_sc <= vmax ||
+                throw(DomainError(si, "sample generated by `[rand(s) for _ in 1:n]` is out of valid range [$vmin, $vmax]."))
         end
     end
 
@@ -174,7 +194,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 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 of samples generated by `[rand(s) for _ in 1:n]` are out of the confidence interval.")
     end
     return samples
 end
@@ -250,13 +274,24 @@ function test_samples(s::Sampleable{Univariate, Continuous},    # the sampleable
         samples = rand(s, n)
         Random.seed!(1234)
         samples2 = rand(s, n)
+        Random.seed!(1234)
+        samples3 = [rand(s) for _ in 1:n]
+        Random.seed!(1234)
+        samples4 = [rand(s) for _ in 1:n]
     else
-        rng2 = 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]
+        samples4 = [rand(rng4, s) for _ in 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.
@@ -266,20 +301,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.")
+            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 ||
+            throw(DomainError(si, "sample generated by `[rand(s) for _ in 1:n]` is out of valid range [$vmin, $vmax]."))
     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 of samples generated by `rand(s, n)` are out of the confidence interval.")
+        clb[i] <= cnts_sc[i] <= cub[i] ||
+            error("The counts of samples generated by `[rand(s) for _ in 1:n]` are out of the confidence interval.")
     end
     return samples
 end
@@ -583,6 +627,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
diff --git a/test/univariate/continuous/exponential.jl b/test/univariate/continuous/exponential.jl
index 191528fd7e..9c5f29aa93 100644
--- a/test/univariate/continuous/exponential.jl
+++ b/test/univariate/continuous/exponential.jl
@@ -1,4 +1,3 @@
-
 @testset "Exponential" begin
     test_cgf(Exponential(1), (0.9, -1, -100f0, -1e6))
     test_cgf(Exponential(0.91), (0.9, -1, -100f0, -1e6))
@@ -8,3 +7,18 @@
         @test @inferred(rand(Exponential(T(1)))) isa T
     end
 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))
+
+# Sampling Tests
+@testset "Exponential sampling tests" begin
+    for d in [
+        Exponential(1),
+        Exponential(0.91),
+        Exponential(10)
+    ]
+        test_distr(d, 10^6)
+    end
+end
diff --git a/test/univariate/continuous/logitnormal.jl b/test/univariate/continuous/logitnormal.jl
index 454376606c..b8d72290ef 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)
+    end
+end
diff --git a/test/univariate/continuous/lognormal.jl b/test/univariate/continuous/lognormal.jl
index 4a27f2e8ac..1cc44b433d 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)
+    end
+end
diff --git a/test/univariate/continuous/normalcanon.jl b/test/univariate/continuous/normalcanon.jl
new file mode 100644
index 0000000000..9ae52e15f0
--- /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)
+    end
+end
diff --git a/test/univariate/continuous/pareto.jl b/test/univariate/continuous/pareto.jl
new file mode 100644
index 0000000000..195e34a96e
--- /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)
+    end
+end

From 61b64b8f3b559d2dda4b6d2f6fbffc362a7dbb20 Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Wed, 25 Sep 2024 13:49:40 +0200
Subject: [PATCH 2/9] Fix deprecation warnings in tests (#1894)

* Fix deprecation warnings in tests

* Change deprecations

* Fixes

* Fix for Julia 1.3

* Fix for Julia 1.3

* `redirect_stderr` on Julia < 1.6
---
 src/Distributions.jl   |   1 -
 src/deprecates.jl      |  15 +++-
 test/matrixreshaped.jl | 189 +++++++++++++++++++++++------------------
 3 files changed, 116 insertions(+), 89 deletions(-)

diff --git a/src/Distributions.jl b/src/Distributions.jl
index 8d344c4158..7cb9f1ed85 100644
--- a/src/Distributions.jl
+++ b/src/Distributions.jl
@@ -129,7 +129,6 @@ export
     MatrixBeta,
     MatrixFDist,
     MatrixNormal,
-    MatrixReshaped,
     MatrixTDist,
     MixtureModel,
     Multinomial,
diff --git a/src/deprecates.jl b/src/deprecates.jl
index 3106154a5e..a360c977c8 100644
--- a/src/deprecates.jl
+++ b/src/deprecates.jl
@@ -53,11 +53,20 @@ end
 @deprecate expectation(distr::Union{UnivariateDistribution,MultivariateDistribution}, g::Function; kwargs...) expectation(g, distr; kwargs...) false
 
 # Deprecate `MatrixReshaped`
+# This is very similar to `Base.@deprecate_binding MatrixReshaped{...} ReshapedDistribution{...}`
+# However, `Base.@deprecate_binding` does not support type parameters
+export MatrixReshaped
 const MatrixReshaped{S<:ValueSupport,D<:MultivariateDistribution{S}} = ReshapedDistribution{2,S,D}
 Base.deprecate(@__MODULE__, :MatrixReshaped)
-@deprecate MatrixReshaped(
-    d::MultivariateDistribution, n::Integer, p::Integer=n
-) reshape(d, (n, p))
+# This is very similar to `Base.@deprecate MatrixReshaped(...) reshape(...)`
+# We use another (unexported!) alias here to not throw a deprecation warning/error
+# Unexported aliases do not affect the type printing
+# In Julia >= 1.6, instead of a new alias we could have defined a method for (ReshapedDistribution{2,S,D} where {S<:ValueSupport,D<:MultivariateDistribution{S}})
+const _MatrixReshaped{S<:ValueSupport,D<:MultivariateDistribution{S}} = ReshapedDistribution{2,S,D}
+function _MatrixReshaped(d::MultivariateDistribution, n::Integer, p::Integer=n)
+    Base.depwarn("`MatrixReshaped(d, n, p)` is deprecated, use `reshape(d, (n, p))` instead.", :MatrixReshaped)
+    return reshape(d, (n, p))
+end
 
 for D in (:InverseWishart, :LKJ, :MatrixBeta, :MatrixFDist, :Wishart)
     @eval @deprecate dim(d::$D) size(d, 1)
diff --git a/test/matrixreshaped.jl b/test/matrixreshaped.jl
index da92ecd2ec..8f5de605fb 100644
--- a/test/matrixreshaped.jl
+++ b/test/matrixreshaped.jl
@@ -3,113 +3,132 @@ using Distributions, Test, Random, LinearAlgebra
 
 rng = MersenneTwister(123456)
 
-@testset "matrixreshaped.jl" begin
+if VERSION >= v"1.6.0-DEV.254"
+    _redirect_stderr(f, ::Base.DevNull) = redirect_stderr(f, devnull)
+else
+    function _redirect_stderr(f, ::Base.DevNull)
+        nulldev = @static Sys.iswindows() ? "NUL" : "/dev/null"
+        open(nulldev, "w") do io
+            redirect_stderr(f, io)
+        end
+    end
+end
 
 function test_matrixreshaped(rng, d1, sizes)
-    x1 = rand(rng, d1)
-    d1s = [@test_deprecated(MatrixReshaped(d1, s...)) for s in sizes]
+    @testset "MatrixReshaped $(nameof(typeof(d1))) tests" begin
+        x1 = rand(rng, d1)
+        d1s = [@test_deprecated(MatrixReshaped(d1, s...)) for s in sizes]
 
-@testset "MatrixReshaped $(nameof(typeof(d1))) tests" begin
-    @testset "MatrixReshaped constructor" begin
-        for d in d1s
-            @test d isa MatrixReshaped
+        @testset "MatrixReshaped constructor" begin
+            for d in d1s
+                @test d isa MatrixReshaped
+            end
         end
-    end
-    @testset "MatrixReshaped constructor errors" begin
-        @test_deprecated(@test_throws ArgumentError MatrixReshaped(d1, length(d1), 2))
-        @test_deprecated(@test_throws ArgumentError MatrixReshaped(d1, length(d1)))
-        @test_deprecated(@test_throws ArgumentError MatrixReshaped(d1, -length(d1), -1))
-    end
-    @testset "MatrixReshaped size" begin
-        for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
-            @test size(d) == s
+        @testset "MatrixReshaped constructor errors" begin
+            @test_deprecated(@test_throws ArgumentError MatrixReshaped(d1, length(d1), 2))
+            @test_deprecated(@test_throws ArgumentError MatrixReshaped(d1, length(d1)))
+            @test_deprecated(@test_throws ArgumentError MatrixReshaped(d1, -length(d1), -1))
         end
-    end
-    @testset "MatrixReshaped length" begin
-        for d in d1s
-            @test length(d) == length(d1)
+        @testset "MatrixReshaped size" begin
+            for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
+                @test size(d) == s
+            end
         end
-    end
-    @testset "MatrixReshaped rank" begin
-        for (d, s) in zip(d1s, sizes)
-            @test rank(d) == minimum(s)
+        @testset "MatrixReshaped length" begin
+            for d in d1s
+                @test length(d) == length(d1)
+            end
         end
-    end
-    @testset "MatrixReshaped insupport" begin
-        for (i, d) in enumerate(d1s[1:end-1])
-            for (j, s) in enumerate(sizes[1:end-1])
-                @test (i == j) ⊻ !insupport(d, reshape(x1, s))
+        @testset "MatrixReshaped rank" begin
+            for (d, s) in zip(d1s, sizes)
+                @test rank(d) == minimum(s)
             end
         end
-    end
-    @testset "MatrixReshaped mean" begin
-        for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
-            @test mean(d) == reshape(mean(d1), s)
+        @testset "MatrixReshaped insupport" begin
+            for (i, d) in enumerate(d1s[1:end-1])
+                for (j, s) in enumerate(sizes[1:end-1])
+                    @test (i == j) ⊻ !insupport(d, reshape(x1, s))
+                end
+            end
         end
-    end
-    @testset "MatrixReshaped mode" begin
-        for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
-            @test mode(d) == reshape(mode(d1), s)
+        @testset "MatrixReshaped mean" begin
+            for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
+                @test mean(d) == reshape(mean(d1), s)
+            end
         end
-    end
-    @testset "MatrixReshaped covariance" begin
-        for (d, (n, p)) in zip(d1s[1:end-1], sizes[1:end-1])
-            @test cov(d) == cov(d1)
-            @test cov(d, Val(false)) == reshape(cov(d1), n, p, n, p)
+        @testset "MatrixReshaped mode" begin
+            for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
+                @test mode(d) == reshape(mode(d1), s)
+            end
         end
-    end
-    @testset "MatrixReshaped variance" begin
-        for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
-            @test var(d) == reshape(var(d1), s)
+        @testset "MatrixReshaped covariance" begin
+            for (d, (n, p)) in zip(d1s[1:end-1], sizes[1:end-1])
+                @test cov(d) == cov(d1)
+                @test cov(d, Val(false)) == reshape(cov(d1), n, p, n, p)
+            end
         end
-    end
-    @testset "MatrixReshaped params" begin
-        for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
-            @test params(d) == (d1, s)
+        @testset "MatrixReshaped variance" begin
+            for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
+                @test var(d) == reshape(var(d1), s)
+            end
         end
-    end
-    @testset "MatrixReshaped partype" begin
-        for d in d1s
-            @test partype(d) === partype(d1)
+        @testset "MatrixReshaped params" begin
+            for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
+                @test params(d) == (d1, s)
+            end
         end
-    end
-    @testset "MatrixReshaped eltype" begin
-        for d in d1s
-            @test eltype(d) === eltype(d1)
+        @testset "MatrixReshaped partype" begin
+            for d in d1s
+                @test partype(d) === partype(d1)
+            end
         end
-    end
-    @testset "MatrixReshaped logpdf" begin
-        for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
-            x = reshape(x1, s)
-            @test logpdf(d, x) == logpdf(d1, x1)
+        @testset "MatrixReshaped eltype" begin
+            for d in d1s
+                @test eltype(d) === eltype(d1)
+            end
         end
-    end
-    @testset "MatrixReshaped rand" begin
-        for d in d1s
-            x = rand(rng, d)
-            @test insupport(d, x)
-            @test insupport(d1, vec(x))
-            @test logpdf(d, x) == logpdf(d1, vec(x))
+        @testset "MatrixReshaped logpdf" begin
+            for (d, s) in zip(d1s[1:end-1], sizes[1:end-1])
+                x = reshape(x1, s)
+                @test logpdf(d, x) == logpdf(d1, x1)
+            end
         end
-    end
-    @testset "MatrixReshaped vec" begin
-        for d in d1s
-            @test vec(d) === d1
+        @testset "MatrixReshaped rand" begin
+            for d in d1s
+                x = rand(rng, d)
+                @test insupport(d, x)
+                @test insupport(d1, vec(x))
+                @test logpdf(d, x) == logpdf(d1, vec(x))
+            end
+        end
+        @testset "MatrixReshaped vec" begin
+            for d in d1s
+                @test vec(d) === d1
+            end
         end
     end
-end
 end
 
-    # MvNormal
-    σ = rand(rng, 16, 16)
-    μ = rand(rng, 16)
-    d1 = MvNormal(μ, σ * σ')
-    sizes = [(4, 4), (8, 2), (2, 8), (1, 16), (16, 1), (4,)]
-    test_matrixreshaped(rng, d1, sizes)
+# Note: In contrast to `@deprecate`, `@deprecate_binding` can't be tested with `@test_deprecated`
+# Ref: https://github.com/JuliaLang/julia/issues/38780
+@testset "matrixreshaped.jl" begin
+    @testset "MvNormal" begin
+        σ = rand(rng, 16, 16)
+        μ = rand(rng, 16)
+        d1 = MvNormal(μ, σ * σ')
+        sizes = [(4, 4), (8, 2), (2, 8), (1, 16), (16, 1), (4,)]
+        _redirect_stderr(devnull) do
+            test_matrixreshaped(rng, d1, sizes)
+        end
+    end
 
     # Dirichlet
-    α = rand(rng, 36) .+ 1 # mode is only defined if all alpha > 1
-    d1 = Dirichlet(α)
-    sizes = [(6, 6), (4, 9), (9, 4), (3, 12), (12, 3), (1, 36), (36, 1), (6,)]
-    test_matrixreshaped(rng, d1, sizes)
+    @testset "Dirichlet" begin
+        α = rand(rng, 36) .+ 1 # mode is only defined if all alpha > 1
+        d1 = Dirichlet(α)
+        sizes = [(6, 6), (4, 9), (9, 4), (3, 12), (12, 3), (1, 36), (36, 1), (6,)]
+        _redirect_stderr(devnull) do
+            test_matrixreshaped(rng, d1, sizes)
+        end
+    end
 end

From a1010e408dbbdd4bcfd9c11eef189df64f7bb05a Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Wed, 25 Sep 2024 13:53:11 +0200
Subject: [PATCH 3/9] Update Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index af12e4b7a1..5cda146036 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Distributions"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
 authors = ["JuliaStats"]
-version = "0.25.111"
+version = "0.25.112"
 
 [deps]
 AliasTables = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"

From aad64af36e83f9a191de34f497e584943ffa84e5 Mon Sep 17 00:00:00 2001
From: Phillip Alday <palday@users.noreply.github.com>
Date: Mon, 4 Nov 2024 21:17:52 +0000
Subject: [PATCH 4/9] Geometric parameter should be on (0, 1] (#1906)

* Geometric parameter should be on (0, 1]

* add tests

* patch bump
---
 Project.toml                          | 2 +-
 src/univariate/discrete/geometric.jl  | 2 +-
 test/univariate/discrete/geometric.jl | 5 +++++
 3 files changed, 7 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 5cda146036..341fb436c1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Distributions"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
 authors = ["JuliaStats"]
-version = "0.25.112"
+version = "0.25.113"
 
 [deps]
 AliasTables = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"
diff --git a/src/univariate/discrete/geometric.jl b/src/univariate/discrete/geometric.jl
index f89cebfedd..f7f940ee4b 100644
--- a/src/univariate/discrete/geometric.jl
+++ b/src/univariate/discrete/geometric.jl
@@ -30,7 +30,7 @@ struct Geometric{T<:Real} <: DiscreteUnivariateDistribution
 end
 
 function Geometric(p::Real; check_args::Bool=true)
-    @check_args Geometric (p, zero(p) < p < one(p))
+    @check_args Geometric (p, zero(p) < p <= one(p))
     return Geometric{typeof(p)}(p)
 end
 
diff --git a/test/univariate/discrete/geometric.jl b/test/univariate/discrete/geometric.jl
index 9845946db7..82f7a1e2ce 100644
--- a/test/univariate/discrete/geometric.jl
+++ b/test/univariate/discrete/geometric.jl
@@ -18,3 +18,8 @@ using FiniteDifferences
     test_cgf(Geometric(0.1), (1f-1, -1e6))
     test_cgf(Geometric(0.5), (1f-1, -1e6))
 end
+
+@testset "Support" begin
+    @test rand(Geometric(1)) == 0
+    @test_throws DomainError Geometric(0)
+end

From 3de6038ef9b041ceaf649bc324d60ca65fcf09eb Mon Sep 17 00:00:00 2001
From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com>
Date: Tue, 19 Nov 2024 09:11:04 +0100
Subject: [PATCH 5/9] Bump codecov/codecov-action from 4 to 5 (#1917)

Bumps [codecov/codecov-action](https://github.com/codecov/codecov-action) from 4 to 5.
- [Release notes](https://github.com/codecov/codecov-action/releases)
- [Changelog](https://github.com/codecov/codecov-action/blob/main/CHANGELOG.md)
- [Commits](https://github.com/codecov/codecov-action/compare/v4...v5)

---
updated-dependencies:
- dependency-name: codecov/codecov-action
  dependency-type: direct:production
  update-type: version-update:semver-major
...

Signed-off-by: dependabot[bot] <support@github.com>
Co-authored-by: dependabot[bot] <49699333+dependabot[bot]@users.noreply.github.com>
---
 .github/workflows/CI.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 886724bd29..2204ab3236 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -47,7 +47,7 @@ jobs:
             Pkg.instantiate()'
       - run: julia --project=perf perf/samplers.jl
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v4
+      - uses: codecov/codecov-action@v5
         with:
           token: ${{ secrets.CODECOV_TOKEN }} # required
           fail_ci_if_error: true

From 2d28655cd9d24f60cd1e5e426e6ad13ff484484b Mon Sep 17 00:00:00 2001
From: Marcus P S <marcusps@users.noreply.github.com>
Date: Sat, 14 Dec 2024 12:19:58 +0000
Subject: [PATCH 6/9] Fix `median` for `Binomial` distribution (#1924)

* Fix `median` for `Binomial`

Fix #1923

* Add tests for `median(x::Binomial)`

* Remove specialization of `median` to `Binomial`

Let `median(x::Binomial)` fall back to `quantile(x,1//2)`

* Update R reference

* Simplify test

---------

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
---
 src/univariate/discrete/binomial.jl  | 2 --
 test/ref/discrete/binomial.R         | 2 +-
 test/ref/discrete_test.ref.json      | 2 +-
 test/univariate/discrete/binomial.jl | 7 +++++++
 4 files changed, 9 insertions(+), 4 deletions(-)

diff --git a/src/univariate/discrete/binomial.jl b/src/univariate/discrete/binomial.jl
index f4102cbb80..5beb4137d0 100644
--- a/src/univariate/discrete/binomial.jl
+++ b/src/univariate/discrete/binomial.jl
@@ -73,8 +73,6 @@ function mode(d::Binomial{T}) where T<:Real
 end
 modes(d::Binomial) = Int[mode(d)]
 
-median(d::Binomial) = round(Int,mean(d))
-
 function skewness(d::Binomial)
     n, p1 = params(d)
     p0 = 1 - p1
diff --git a/test/ref/discrete/binomial.R b/test/ref/discrete/binomial.R
index 31d85a0064..b6d8e14fac 100644
--- a/test/ref/discrete/binomial.R
+++ b/test/ref/discrete/binomial.R
@@ -18,7 +18,7 @@ Binomial <- R6Class("Binomial",
                  failprob=q,
                  ntrials=n,
                  mean=n * p,
-                 median=round(n * p),
+                 median=qbinom(0.5, self$n, self$p),
                  var=n * p * q,
                  skewness=(q - p) / sqrt(n*p*q),
                  kurtosis=(1 - 6*p*q) / (n*p*q))
diff --git a/test/ref/discrete_test.ref.json b/test/ref/discrete_test.ref.json
index 17f9cd2d51..05a465be71 100644
--- a/test/ref/discrete_test.ref.json
+++ b/test/ref/discrete_test.ref.json
@@ -285,7 +285,7 @@
     "failprob": 0.5,
     "ntrials": 3,
     "mean": 1.5,
-    "median": 2,
+    "median": 1,
     "var": 0.75,
     "skewness": 0,
     "kurtosis": -0.666666666666667
diff --git a/test/univariate/discrete/binomial.jl b/test/univariate/discrete/binomial.jl
index 0caa974216..5f95735bbb 100644
--- a/test/univariate/discrete/binomial.jl
+++ b/test/univariate/discrete/binomial.jl
@@ -26,6 +26,13 @@ end
 @test Distributions.expectation(identity, Binomial(6)) ≈ 3.0
 @test Distributions.expectation(x -> -x, Binomial(10, 0.2)) ≈ -2.0
 
+# Test median
+@test median(Binomial(5,3//10)) == 1
+@test median(Binomial(25,3//10)) == 7
+@test median(Binomial(45,3//10)) == 13
+@test median(Binomial(65,3//10)) == 19
+@test median(Binomial(85,3//10)) == 25
+    
 # Test mode
 @test Distributions.mode(Binomial(100, 0.4)) == 40
 @test Distributions.mode(Binomial(1, 0.51)) == 1

From 790411a14f66311bd5caada76105958d4fa542a8 Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Sat, 14 Dec 2024 13:20:25 +0100
Subject: [PATCH 7/9] Update Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 341fb436c1..5c0f8dd15c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Distributions"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
 authors = ["JuliaStats"]
-version = "0.25.113"
+version = "0.25.114"
 
 [deps]
 AliasTables = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"

From 445c2fe317c90533a08d8ba69d17c31b336a75f9 Mon Sep 17 00:00:00 2001
From: Marcus P S <marcusps@users.noreply.github.com>
Date: Tue, 17 Dec 2024 11:17:01 +0000
Subject: [PATCH 8/9] Re-write median(Binomial) to improve speed (#1928)

* Speed up median(Binomial)

Follow suggestion from #1926 to reduce `cdf(Binomial)` evaluations.

* Simplify bound condition

Co-authored-by: Andreas Noack <andreas@noack.dk>

* Fix name of argument variable

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>

* Simplify b

Co-authored-by: Andreas Noack <andreas@noack.dk>

* Fix cummulative probability check

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>

* Incorporate suggestions from PR review

* Further simplifications following suggestions

* Add more tests

* Fixed checks for median bounds

* Simplified tests a bit

* Update src/univariate/discrete/binomial.jl

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>

* Update src/univariate/discrete/binomial.jl

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>

* Update src/univariate/discrete/binomial.jl

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>

* Avoid test memory allocation with an iterator

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>

---------

Co-authored-by: Andreas Noack <andreas@noack.dk>
Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
---
 src/univariate/discrete/binomial.jl  | 31 ++++++++++++++++++++++++++++
 test/univariate/discrete/binomial.jl |  3 ++-
 2 files changed, 33 insertions(+), 1 deletion(-)

diff --git a/src/univariate/discrete/binomial.jl b/src/univariate/discrete/binomial.jl
index 5beb4137d0..74c90ee18c 100644
--- a/src/univariate/discrete/binomial.jl
+++ b/src/univariate/discrete/binomial.jl
@@ -73,6 +73,37 @@ function mode(d::Binomial{T}) where T<:Real
 end
 modes(d::Binomial) = Int[mode(d)]
 
+function median(dist::Binomial)
+    # The median is floor(Int, mean) or ceil(Int, mean)
+    # As shown in https://doi.org/10.1016/0167-7152(94)00090-U,
+    # |median - mean| <= 1 - bound
+    # where the equality is strict except for the case p = 1/2 and n odd.
+    # Thus if |k - mean| < bound for one of the two candidates if p = 1/2 and n odd
+    # or |k - mean| <= bound for one of the two candidates otherwise,
+    # the other candidate can't satisfy the condition and hence k must be the median
+    bound = max(min(dist.p, 1-dist.p), loghalf)
+    dist_mean = mean(dist)
+    
+    floor_mean = floor(Int, dist_mean)
+    difference = dist_mean - floor_mean
+    
+    if difference <= bound
+        # The only case where the median satisfies |median - mean| <= 1 - bound with equality
+        # is p = 1/2 and n odd
+        # However, in that case we also want to return floor(mean)
+        floor_mean
+    elseif difference >= 1 - bound
+        # The case p = 1/2 and n odd was already covered above,
+        # thus only cases with |median - mean| < 1 - bound are left here
+        # Therefore difference >= 1 - bound implies that floor(mean) cannot be the median
+        floor_mean + 1
+    elseif cdf(dist, floor_mean) >= 0.5
+        floor_mean
+    else
+        floor_mean + 1
+    end
+end
+
 function skewness(d::Binomial)
     n, p1 = params(d)
     p0 = 1 - p1
diff --git a/test/univariate/discrete/binomial.jl b/test/univariate/discrete/binomial.jl
index 5f95735bbb..55b3362144 100644
--- a/test/univariate/discrete/binomial.jl
+++ b/test/univariate/discrete/binomial.jl
@@ -27,11 +27,12 @@ end
 @test Distributions.expectation(x -> -x, Binomial(10, 0.2)) ≈ -2.0
 
 # Test median
-@test median(Binomial(5,3//10)) == 1
 @test median(Binomial(25,3//10)) == 7
 @test median(Binomial(45,3//10)) == 13
 @test median(Binomial(65,3//10)) == 19
 @test median(Binomial(85,3//10)) == 25
+
+@test all(median(Binomial(7, p)) == quantile(Binomial(7, p), 1//2) for p in 0:0.1:1)
     
 # Test mode
 @test Distributions.mode(Binomial(100, 0.4)) == 40

From ceb63433f0624881484b499cf05e664b7a446b07 Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Tue, 17 Dec 2024 12:17:21 +0100
Subject: [PATCH 9/9] Update Project.toml

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 5c0f8dd15c..ac4358fb0b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "Distributions"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
 authors = ["JuliaStats"]
-version = "0.25.114"
+version = "0.25.115"
 
 [deps]
 AliasTables = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"