From b09dcd797ec5f191592194718e3b34ff55cb7963 Mon Sep 17 00:00:00 2001
From: Carlos Parada <paradac@carleton.edu>
Date: Fri, 17 May 2024 14:21:52 -0700
Subject: [PATCH 01/19] Fix `fit` docs (#1746)

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
Co-authored-by: Viral B. Shah <ViralBShah@users.noreply.github.com>
---
 docs/src/fit.md                   |  2 ++
 src/genericfit.jl                 | 13 +++++++++++++
 src/univariate/continuous/beta.jl | 11 ++---------
 3 files changed, 17 insertions(+), 9 deletions(-)

diff --git a/docs/src/fit.md b/docs/src/fit.md
index 24ddbad5ce..b482f995e6 100644
--- a/docs/src/fit.md
+++ b/docs/src/fit.md
@@ -27,6 +27,8 @@ The function `fit_mle` is for maximum likelihood estimation.
 ### Synopsis
 
 ```@docs
+fit(D, x)
+fit(D, x, w)
 fit_mle(D, x)
 fit_mle(D, x, w)
 ```
diff --git a/src/genericfit.jl b/src/genericfit.jl
index 2b4b5e67eb..b6bce07a6a 100644
--- a/src/genericfit.jl
+++ b/src/genericfit.jl
@@ -30,5 +30,18 @@ fit_mle(dt::Type{D}, x::AbstractArray, w::AbstractArray) where {D<:UnivariateDis
 fit_mle(dt::Type{D}, x::AbstractMatrix) where {D<:MultivariateDistribution} = fit_mle(D, suffstats(D, x))
 fit_mle(dt::Type{D}, x::AbstractMatrix, w::AbstractArray) where {D<:MultivariateDistribution} = fit_mle(D, suffstats(D, x, w))
 
+"""
+    fit(D, args...)
+
+Fit a distribution of type `D` to `args`.
+
+The fit function will choose a reasonable way to fit the distribution, which,
+in most cases, is maximum likelihood estimation. Note that this algorithm may
+change; for a function that will behave consistently across versions, see 
+`fit_mle`.
+
+By default, the fallback is [`fit_mle(D, args...)`](@ref); developers can change this default
+for a specific distribution type `D <: Distribution` by defining a `fit(::Type{D}, args...)` method.
+"""
 fit(dt::Type{D}, x) where {D<:Distribution} = fit_mle(D, x)
 fit(dt::Type{D}, args...) where {D<:Distribution} = fit_mle(D, args...)
diff --git a/src/univariate/continuous/beta.jl b/src/univariate/continuous/beta.jl
index 44501fb3c3..fd2420d815 100644
--- a/src/univariate/continuous/beta.jl
+++ b/src/univariate/continuous/beta.jl
@@ -208,12 +208,8 @@ function rand(rng::AbstractRNG, d::Beta{T}) where T
     end
 end
 
-#### Fit model
-"""
-    fit_mle(::Type{<:Beta}, x::AbstractArray{T})
 
-Maximum Likelihood Estimate of `Beta` Distribution via Newton's Method
-"""
+
 function fit_mle(::Type{<:Beta}, x::AbstractArray{T};
     maxiter::Int=1000, tol::Float64=1e-14) where T<:Real
 
@@ -240,11 +236,8 @@ function fit_mle(::Type{<:Beta}, x::AbstractArray{T};
     return Beta(θ[1], θ[2])
 end
 
-"""
-    fit(::Type{<:Beta}, x::AbstractArray{T})
 
-fit a `Beta` distribution
-"""
+
 function fit(::Type{<:Beta}, x::AbstractArray{T}) where T<:Real
     x_bar = mean(x)
     v_bar = varm(x, x_bar)

From 04c7bf61fdb26eae556d5422b09bf49e18980331 Mon Sep 17 00:00:00 2001
From: "Viral B. Shah" <ViralBShah@users.noreply.github.com>
Date: Sun, 19 May 2024 00:16:24 -0400
Subject: [PATCH 02/19] Update CI.yml to include mac aarch64

---
 .github/workflows/CI.yml | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index ccf2aea36f..9d791e2056 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -33,6 +33,10 @@ jobs:
           - windows-latest
         arch:
           - x64
+        include:
+          - os: macOS-14
+            arch: aarch64
+            version: '1'
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2

From 24de00402125d876ec378ceeb6fbe408bfecf75a Mon Sep 17 00:00:00 2001
From: Tim Hargreaves <38204689+THargreaves@users.noreply.github.com>
Date: Wed, 29 May 2024 21:51:04 +0100
Subject: [PATCH 03/19] Replace depricated lfact(x) (#1857)

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

diff --git a/src/univariate/continuous/ksdist.jl b/src/univariate/continuous/ksdist.jl
index 674bcec95a..1de09651a3 100644
--- a/src/univariate/continuous/ksdist.jl
+++ b/src/univariate/continuous/ksdist.jl
@@ -28,7 +28,7 @@ function cdf(d::KSDist,x::Float64)
         return 0.0
     elseif b <= 1
         # accuracy could be improved
-        return exp(lfact(n)+n*(log(2*b-1)-log(n)))
+        return exp(logfactorial(n)+n*(log(2*b-1)-log(n)))
     elseif x >= 1
         return 1.0
     elseif b >= n-1
@@ -56,7 +56,7 @@ function ccdf(d::KSDist,x::Float64)
     if b <= 0.5
         return 1.0
     elseif b <= 1
-        return 1-exp(lfact(n)+n*(log(2*b-1)-log(n)))
+        return 1-exp(logfactorial(n)+n*(log(2*b-1)-log(n)))
     elseif x >= 1
         return 0.0
     elseif b >= n-1

From fe57164dbbe710ead723b22e3cd01c9f0870311a Mon Sep 17 00:00:00 2001
From: ajshephard <andrew.james.shephard@gmail.com>
Date: Wed, 29 May 2024 16:53:07 -0400
Subject: [PATCH 04/19] fixing error in xval function in gumbel.jl (#1859)

* fixing error in xval function in gumbel.jl

* implement Gumbel quantile function using xval

* adding Gumbel functions: ccdf, logccdf, cquantile, invlogcdf invlogccdf, mgf, cgf, cf

* Update src/univariate/continuous/gumbel.jl

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

* Update src/univariate/continuous/gumbel.jl

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

* Update src/univariate/continuous/gumbel.jl

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

---------

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
---
 src/univariate/continuous/gumbel.jl | 13 +++++++++++--
 1 file changed, 11 insertions(+), 2 deletions(-)

diff --git a/src/univariate/continuous/gumbel.jl b/src/univariate/continuous/gumbel.jl
index 7977a17027..54090967fc 100644
--- a/src/univariate/continuous/gumbel.jl
+++ b/src/univariate/continuous/gumbel.jl
@@ -85,7 +85,7 @@ entropy(d::Gumbel) = log(d.θ) + 1 + MathConstants.γ
 #### Evaluation
 
 zval(d::Gumbel, x::Real) = (x - d.μ) / d.θ
-xval(d::Gumbel, z::Real) = x * d.θ + d.μ
+xval(d::Gumbel, z::Real) = z * d.θ + d.μ
 
 function pdf(d::Gumbel, x::Real)
     z = zval(d, x)
@@ -98,8 +98,17 @@ function logpdf(d::Gumbel, x::Real)
 end
 
 cdf(d::Gumbel, x::Real) = exp(-exp(-zval(d, x)))
+ccdf(d::Gumbel, x::Real) = -expm1(-exp(-zval(d, x)))
 logcdf(d::Gumbel, x::Real) = -exp(-zval(d, x))
+logccdf(d::Gumbel, x::Real) = log1mexp(-exp(-zval(d, x)))
 
-quantile(d::Gumbel, p::Real) = d.μ - d.θ * log(-log(p))
+quantile(d::Gumbel, p::Real) = xval(d, -log(-log(p)))
+cquantile(d::Gumbel, p::Real) = xval(d, -log(-log1p(-p)))
+invlogcdf(d::Gumbel, lp::Real) = xval(d, -log(-lp))
+invlogccdf(d::Gumbel, lp::Real) = xval(d, -log(-log1mexp(lp)))
 
 gradlogpdf(d::Gumbel, x::Real) = expm1(-zval(d, x)) / d.θ
+
+mgf(d::Gumbel, t::Real) = gamma(1 - d.θ * t) * exp(d.μ * t)
+cgf(d::Gumbel, t::Real) = loggamma(1 - d.θ * t) + d.μ * t
+cf(d::Gumbel, t::Real) = gamma(1 - im * d.θ * t) * cis(d.μ * t)

From fa493fb84020b9c8f9449ffeccbc5b0d6c313831 Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Thu, 30 May 2024 08:58:28 +0200
Subject: [PATCH 05/19] Clamp `cdf` and `ccdf` of `Truncated` (#1865)

---
 src/truncate.jl  |  6 ++++--
 test/truncate.jl | 18 ++++++++++++++++++
 2 files changed, 22 insertions(+), 2 deletions(-)

diff --git a/src/truncate.jl b/src/truncate.jl
index 45709f6b53..48d62b015b 100644
--- a/src/truncate.jl
+++ b/src/truncate.jl
@@ -166,7 +166,8 @@ function logpdf(d::Truncated, x::Real)
 end
 
 function cdf(d::Truncated, x::Real)
-    result = (cdf(d.untruncated, x) - d.lcdf) / d.tp
+    result = clamp((cdf(d.untruncated, x) - d.lcdf) / d.tp, 0, 1)
+    # Special cases for values outside of the support to avoid e.g. NaN issues with `Binomial`
     return if d.lower !== nothing && x < d.lower
         zero(result)
     elseif d.upper !== nothing && x >= d.upper
@@ -188,7 +189,8 @@ function logcdf(d::Truncated, x::Real)
 end
 
 function ccdf(d::Truncated, x::Real)
-    result = (d.ucdf - cdf(d.untruncated, x)) / d.tp
+    result = clamp((d.ucdf - cdf(d.untruncated, x)) / d.tp, 0, 1)
+    # Special cases for values outside of the support to avoid e.g. NaN issues with `Binomial`
     return if d.lower !== nothing && x <= d.lower
         one(result)
     elseif d.upper !== nothing && x > d.upper
diff --git a/test/truncate.jl b/test/truncate.jl
index b409cba93e..9c2a286d09 100644
--- a/test/truncate.jl
+++ b/test/truncate.jl
@@ -214,3 +214,21 @@ end
 
   @test isa(quantile(d, ForwardDiff.Dual(1.,0.)), ForwardDiff.Dual)
 end
+
+@testset "cdf outside of [0, 1] (#1854)" begin
+    dist = truncated(Normal(2.5, 0.2); lower=0.0)
+    @test @inferred(cdf(dist, 3.741058503233821e-17)) === 0.0
+    @test @inferred(ccdf(dist, 3.741058503233821e-17)) === 1.0
+    @test @inferred(cdf(dist, 1.4354474178676617e-18)) === 0.0
+    @test @inferred(ccdf(dist, 1.4354474178676617e-18)) === 1.0
+    @test @inferred(cdf(dist, 8.834854780587132e-18)) === 0.0
+    @test @inferred(ccdf(dist, 8.834854780587132e-18)) === 1.0
+
+    dist = truncated(
+        Normal(2.122039143928797, 0.07327367710864985);
+        lower = 1.9521656132878236,
+        upper = 2.8274429454898398,
+    )
+    @test @inferred(cdf(dist, 2.82)) === 1.0
+    @test @inferred(ccdf(dist, 2.82)) === 0.0
+end

From 6af1e2f2f2448006e308f4ab95ca0f66037563e6 Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Thu, 30 May 2024 09:52:42 +0200
Subject: [PATCH 06/19] Update Project.toml

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

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

From b356da03a189d023cdb8467c61806a8a11dcb262 Mon Sep 17 00:00:00 2001
From: Michael Harradon <mharradon@gmail.com>
Date: Thu, 30 May 2024 04:28:22 -0400
Subject: [PATCH 07/19] Add Convolve for DiscreteNonParametric (Redux) (#1850)

* Add convolve for DiscreteNonParametric

DiscreteNonParametric convolution has a very nice trivial closed form. It was not implemented.

This pull request implements it.

* Update src/convolution.jl

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

* Use Set, instead of splatting.

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

* Fix type stability of elements.

Doesn't preserve the type of the Vector, but perhaps this is better ....

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

* Apply suggestions from code review

use functions to access the support and probabilities, and write as one loop.

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

* Added a test set.
Removed check args:
We know the convovultion is a proper distribution.

* minor rename for consistency

* Formatting, test improvements suggested by devmotion (and a few more)

* Formatting

---------

Co-authored-by: iampritishpatil <iampritishpatil@gmail.com>
Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
---
 src/convolution.jl  | 14 ++++++++++++++
 test/convolution.jl | 22 ++++++++++++++++++++++
 2 files changed, 36 insertions(+)

diff --git a/src/convolution.jl b/src/convolution.jl
index 6c18b70e48..cd429e8b42 100644
--- a/src/convolution.jl
+++ b/src/convolution.jl
@@ -12,6 +12,7 @@ and one of
 * [`NegativeBinomial`](@ref)
 * [`Geometric`](@ref)
 * [`Poisson`](@ref)
+* [`DiscreteNonParametric`](@ref)
 * [`Normal`](@ref)
 * [`Cauchy`](@ref)
 * [`Chisq`](@ref)
@@ -47,6 +48,19 @@ end
 convolve(d1::Poisson, d2::Poisson) =  Poisson(d1.λ + d2.λ)
 
 
+function convolve(d1::DiscreteNonParametric, d2::DiscreteNonParametric)
+    support_conv = collect(Set(s1 + s2 for s1 in support(d1), s2 in support(d2)))
+    sort!(support_conv) #for fast index finding below
+    probs1 = probs(d1)
+    probs2 = probs(d2)
+    p_conv = zeros(Base.promote_eltype(probs1, probs2), length(support_conv)) 
+    for (s1, p1) in zip(support(d1), probs(d1)), (s2, p2) in zip(support(d2), probs(d2))
+            idx = searchsortedfirst(support_conv, s1+s2)
+            p_conv[idx] += p1*p2
+    end
+    DiscreteNonParametric(support_conv, p_conv,check_args=false) 
+end
+
 # continuous univariate
 convolve(d1::Normal, d2::Normal) = Normal(d1.μ + d2.μ, hypot(d1.σ, d2.σ))
 convolve(d1::Cauchy, d2::Cauchy) = Cauchy(d1.μ + d2.μ, d1.σ + d2.σ)
diff --git a/test/convolution.jl b/test/convolution.jl
index 826e4f82ce..38ea3943d4 100644
--- a/test/convolution.jl
+++ b/test/convolution.jl
@@ -66,6 +66,28 @@ using Test
         @test d3 isa Poisson
         @test d3.λ == 0.5
     end
+
+    @testset "DiscreteNonParametric" begin
+        d1 = DiscreteNonParametric([0, 1], [0.5, 0.5])
+        d2 = DiscreteNonParametric([1, 2], [0.5, 0.5])
+        d_eps = DiscreteNonParametric([prevfloat(0.0), 0.0, nextfloat(0.0), 1.0], fill(1//4, 4))
+        d10 = DiscreteNonParametric((1//10):(1//10):1, fill(1//10, 10))
+        
+        d_int_simple = @inferred(convolve(d1, d2))
+        @test d_int_simple isa DiscreteNonParametric
+        @test support(d_int_simple) == [1, 2, 3]
+        @test probs(d_int_simple) == [0.25, 0.5, 0.25]
+               
+        d_rat = convolve(d10, d10)
+        @test support(d_rat) == (1//5):(1//10):2
+        @test probs(d_rat) == [1//100, 1//50, 3//100, 1//25, 1//20, 3//50, 7//100, 2//25, 9//100, 1//10, 
+                               9//100, 2//25, 7//100, 3//50, 1//20, 1//25, 3//100, 1//50, 1//100]
+
+        d_float_supp = convolve(d_eps, d_eps)
+        @test support(d_float_supp) == [2 * prevfloat(0.0), prevfloat(0.0), 0.0, nextfloat(0.0), 2 * nextfloat(0.0), 1.0, 2.0]
+        @test probs(d_float_supp) == [1//16, 1//8, 3//16, 1//8, 1//16, 3//8, 1//16]
+    end
+    
 end
 
 @testset "continuous univariate" begin

From 65f056c1f2db29e9a372cdd66b3490b633c4f7e9 Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Sat, 29 Jun 2024 00:06:58 +0200
Subject: [PATCH 08/19] Update CI: No nightly tests, improve macOS CI setup +
 use julia-actions/cache (#1863)

* Test prereleases if available

* Fix name?

* Fix matrix

* Fixes

* Fix syntax

* Another round of fixes

* Remove prerelease CI

* Better mac fix?

* Fix typo

* Fix failing tests on 1.11-beta (#1862)

* Increase sample size when testing empirical moments of DiscreteNonParametric

* Use StableRNGs for semicircle.jl to avoid breakage on 1.11-beta

* Run CI on all PRs. Also CI for workflow_dispatch and merge_group

* Update semicircle.jl

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

* Update locationscale.jl

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

* Update locationscale.jl

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

---------

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

---------

Co-authored-by: Andreas Noack <andreas@noack.dk>
---
 .github/workflows/CI.yml                 | 30 +++++++-----------------
 test/univariate/continuous/semicircle.jl |  6 ++---
 test/univariate/locationscale.jl         | 14 +++++------
 3 files changed, 18 insertions(+), 32 deletions(-)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 9d791e2056..3039dbfc98 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -2,12 +2,12 @@ name: CI
 
 on:
   pull_request:
-    branches:
-      - master
   push:
     branches:
       - master
     tags: '*'
+  workflow_dispatch:
+  merge_group:
 
 concurrency:
   # Skip intermediate builds: always.
@@ -17,42 +17,27 @@ concurrency:
 
 jobs:
   test:
-    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ github.event_name }}
     runs-on: ${{ matrix.os }}
-    continue-on-error: ${{ matrix.version == 'nightly' }}
     strategy:
       fail-fast: false
       matrix:
         version:
           - '1.3'
           - '1'
-          - 'nightly'
         os:
           - ubuntu-latest
           - macos-latest
           - windows-latest
-        arch:
-          - x64
-        include:
-          - os: macOS-14
-            arch: aarch64
-            version: '1'
     steps:
       - uses: actions/checkout@v4
       - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
-          arch: ${{ matrix.arch }}
-      - uses: actions/cache@v4
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-
-            ${{ runner.os }}-
+          # ARM64 on macos-latest is neither supported by older Julia versions nor setup-julia
+          arch: ${{ matrix.os == 'macos-latest' && matrix.version != '1.3' && 'aarch64' || 'x64' }}
+          show-versioninfo: true
+      - uses: julia-actions/cache@v1
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - run: |
@@ -74,6 +59,7 @@ jobs:
       - uses: julia-actions/setup-julia@v2
         with:
           version: '1'
+          show-versioninfo: true
       - run: |
           julia --project=docs -e '
             using Pkg
diff --git a/test/univariate/continuous/semicircle.jl b/test/univariate/continuous/semicircle.jl
index 8dbf2ff2c4..7079aa9490 100644
--- a/test/univariate/continuous/semicircle.jl
+++ b/test/univariate/continuous/semicircle.jl
@@ -1,5 +1,5 @@
 using Distributions
-using Random: MersenneTwister
+using StableRNGs
 using Test
 
 d = Semicircle(2.0)
@@ -39,8 +39,8 @@ d = Semicircle(2.0)
 @test quantile(d,  .5) ==   .0
 @test quantile(d, 1.0) == +2.0
 
-rng = MersenneTwister(0)
-for r in rand(rng, Uniform(0,10), 5)
+rng = StableRNG(123)
+@testset for r in rand(rng, Uniform(0,10), 5)
     N = 10^4
     semi = Semicircle(r)
     sample = rand(rng, semi, N)
diff --git a/test/univariate/locationscale.jl b/test/univariate/locationscale.jl
index c4876e1970..5d5f9640b4 100644
--- a/test/univariate/locationscale.jl
+++ b/test/univariate/locationscale.jl
@@ -110,7 +110,7 @@ function test_location_scale(
             @test invlogccdf(dtest, log(0.5)) ≈ invlogccdf(dref, log(0.5))
             @test invlogccdf(dtest, log(0.8)) ≈ invlogccdf(dref, log(0.8))
 
-            r = Array{float(eltype(dtest))}(undef, 100000)
+            r = Array{float(eltype(dtest))}(undef, 200000)
             if ismissing(rng)
                 rand!(dtest, r)
             else
@@ -148,7 +148,7 @@ end
     rng = MersenneTwister(123)
 
     @testset "Normal" begin
-        for _rng in (missing, rng), sign in (1, -1)
+        @testset for _rng in (missing, rng), sign in (1, -1)
             test_location_scale_normal(_rng, 0.3, sign * 0.2, 0.1, 0.2)
             test_location_scale_normal(_rng, -0.3, sign * 0.1, -0.1, 0.3)
             test_location_scale_normal(_rng, 1.3, sign * 0.4, -0.1, 0.5)
@@ -156,11 +156,11 @@ end
         test_location_scale_normal(rng, ForwardDiff.Dual(0.3), 0.2, 0.1, 0.2)
     end
     @testset "DiscreteNonParametric" begin
-        probs = normalize!(rand(10), 1)
-        for _rng in (missing, rng), sign in (1, -1)
-            test_location_scale_discretenonparametric(_rng, 1//3, sign * 1//2, 1:10, probs)
-            test_location_scale_discretenonparametric(_rng, -1//4, sign * 1//3, (-10):(-1), probs)
-            test_location_scale_discretenonparametric(_rng, 6//5, sign * 3//2, 15:24, probs)
+        _probs = normalize!(rand(10), 1)
+        @testset for _rng in (missing, rng), sign in (1, -1)
+            test_location_scale_discretenonparametric(_rng, 1//3, sign * 1//2, 1:10, _probs)
+            test_location_scale_discretenonparametric(_rng, -1//4, sign * 1//3, (-10):(-1), _probs)
+            test_location_scale_discretenonparametric(_rng, 6//5, sign * 3//2, 15:24, _probs)
         end
     end
 

From 47c040beef8c61bad3e1eefa4fc8194e3a62b55a Mon Sep 17 00:00:00 2001
From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com>
Date: Mon, 1 Jul 2024 21:10:10 +0200
Subject: [PATCH 09/19] Bump julia-actions/cache from 1 to 2 (#1876)

Bumps [julia-actions/cache](https://github.com/julia-actions/cache) from 1 to 2.
- [Release notes](https://github.com/julia-actions/cache/releases)
- [Changelog](https://github.com/julia-actions/cache/blob/main/devdocs/making_a_new_release.md)
- [Commits](https://github.com/julia-actions/cache/compare/v1...v2)

---
updated-dependencies:
- dependency-name: julia-actions/cache
  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 3039dbfc98..384a02a315 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -37,7 +37,7 @@ jobs:
           # ARM64 on macos-latest is neither supported by older Julia versions nor setup-julia
           arch: ${{ matrix.os == 'macos-latest' && matrix.version != '1.3' && 'aarch64' || 'x64' }}
           show-versioninfo: true
-      - uses: julia-actions/cache@v1
+      - uses: julia-actions/cache@v2
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - run: |

From 6e0f1dcec087b517db8f1b62276a3027ddcce5d2 Mon Sep 17 00:00:00 2001
From: quildtide <42811940+quildtide@users.noreply.github.com>
Date: Fri, 9 Aug 2024 05:38:02 -0400
Subject: [PATCH 10/19] Fix InverseGaussian cdf overflows (#1882)

* Fix inversegaussian cdf overflows

* Update src/univariate/continuous/inversegaussian.jl

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

* Update src/univariate/continuous/inversegaussian.jl

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

* Fix missing parentheses

* Move parentheses

* Move new tests to dedicated inversegaussian file

* Add ccdf tests

---------

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
---
 src/univariate/continuous/inversegaussian.jl  | 10 ++++++++--
 test/runtests.jl                              |  2 +-
 test/univariate/continuous/inversegaussian.jl | 18 ++++++++++++++++++
 3 files changed, 27 insertions(+), 3 deletions(-)
 create mode 100644 test/univariate/continuous/inversegaussian.jl

diff --git a/src/univariate/continuous/inversegaussian.jl b/src/univariate/continuous/inversegaussian.jl
index a076d85b21..b585ec3b3d 100644
--- a/src/univariate/continuous/inversegaussian.jl
+++ b/src/univariate/continuous/inversegaussian.jl
@@ -99,7 +99,10 @@ function cdf(d::InverseGaussian, x::Real)
     y = max(x, 0)
     u = sqrt(λ / y)
     v = y / μ
-    z = normcdf(u * (v - 1)) + exp(2λ / μ) * normcdf(-u * (v + 1))
+    # 2λ/μ and normlogcdf(-u*(v+1)) are similar magnitude, opp. sign
+    # truncating to [0, 1] as an additional precaution
+    # Ref https://github.com/JuliaStats/Distributions.jl/issues/1873
+    z = clamp(normcdf(u * (v - 1)) + exp(2λ / μ + normlogcdf(-u * (v + 1))), 0, 1)
 
     # otherwise `NaN` is returned for `+Inf`
     return isinf(x) && x > 0 ? one(z) : z
@@ -110,7 +113,10 @@ function ccdf(d::InverseGaussian, x::Real)
     y = max(x, 0)
     u = sqrt(λ / y)
     v = y / μ
-    z = normccdf(u * (v - 1)) - exp(2λ / μ) * normcdf(-u * (v + 1))
+    # 2λ/μ and normlogcdf(-u*(v+1)) are similar magnitude, opp. sign
+    # truncating to [0, 1] as an additional precaution
+    # Ref https://github.com/JuliaStats/Distributions.jl/issues/1873
+    z = clamp(normccdf(u * (v - 1)) - exp(2λ / μ + normlogcdf(-u * (v + 1))), 0, 1)
 
     # otherwise `NaN` is returned for `+Inf`
     return isinf(x) && x > 0 ? zero(z) : z
diff --git a/test/runtests.jl b/test/runtests.jl
index 583132c536..c7d73f4546 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -96,6 +96,7 @@ const tests = [
     "eachvariate",
     "univariate/continuous/triangular",
     "statsapi",
+    "univariate/continuous/inversegaussian",
 
     ### missing files compared to /src:
     # "common",
@@ -137,7 +138,6 @@ const tests = [
     # "univariate/continuous/generalizedextremevalue",
     # "univariate/continuous/generalizedpareto",
     # "univariate/continuous/inversegamma",
-    # "univariate/continuous/inversegaussian",
     # "univariate/continuous/ksdist",
     # "univariate/continuous/ksonesided",
     # "univariate/continuous/levy",
diff --git a/test/univariate/continuous/inversegaussian.jl b/test/univariate/continuous/inversegaussian.jl
new file mode 100644
index 0000000000..6990cae710
--- /dev/null
+++ b/test/univariate/continuous/inversegaussian.jl
@@ -0,0 +1,18 @@
+@testset "InverseGaussian cdf outside of [0, 1] (#1873)" begin
+    for d in [
+        InverseGaussian(1.65, 590),
+        InverseGaussian(0.5, 1000)
+    ]
+        for x in [0.02, 1.0, 20.0, 300.0]
+            p = cdf(d, x)
+            @test 0.0 <= p <= 1.0
+            @test p ≈ exp(logcdf(d, x))
+
+            q = ccdf(d, x)
+            @test 0.0 <= q <= 1.0
+            @test q ≈ exp(logccdf(d, x))
+
+            @test (p + q) ≈ 1
+        end
+    end
+end
\ No newline at end of file

From 13029c03fa885a2b4512b954e61f9d5a7dfa0612 Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Fri, 9 Aug 2024 11:38:22 +0200
Subject: [PATCH 11/19] Update Project.toml

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

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

From 3946acc6901b5126022ae22564ab9aa7d913ed1e Mon Sep 17 00:00:00 2001
From: Kyurae Kim <msca8h@naver.com>
Date: Fri, 23 Aug 2024 15:44:22 -0400
Subject: [PATCH 12/19] fix type stability of sampling from `Chisq`, `TDist`,
 `Gamma` (#1885)

* fix type stability of sampling from `Chisq`, `TDist`, `Gamma`

* fix remove type specification in `rand(Exponential)`

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

* fix type specificaton in `rand(TDist)`

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

* fix remove type test for `rand(Chisq)`

* fix make `Exponential` use the `Normal` sampling type policy

* fix missing type signature

* fix type signature for `rand(Exponential)`

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

* fix use `@inferred` in tests for `Gamma`

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

* fix use `@inferred` in tests for `TDist`

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

* add type stability tests for `rand(Exponential)`

* add type stability test for `rand(Chisq)`

* fix remove type stability test for `entropy(TDist)` (not stable)

---------

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
---
 src/samplers/gamma.jl                     |  2 +-
 src/univariate/continuous/exponential.jl  |  2 +-
 src/univariate/continuous/tdist.jl        |  2 +-
 test/univariate/continuous/chisq.jl       | 11 +++++--
 test/univariate/continuous/exponential.jl | 12 +++++--
 test/univariate/continuous/gamma.jl       | 38 ++++++++++++++---------
 test/univariate/continuous/tdist.jl       | 17 +++++++---
 7 files changed, 56 insertions(+), 28 deletions(-)

diff --git a/src/samplers/gamma.jl b/src/samplers/gamma.jl
index 9a40da98a3..1cd62f22ba 100644
--- a/src/samplers/gamma.jl
+++ b/src/samplers/gamma.jl
@@ -225,6 +225,6 @@ end
 
 function rand(rng::AbstractRNG, s::GammaIPSampler)
     x = rand(rng, s.s)
-    e = randexp(rng)
+    e = randexp(rng, typeof(x))
     x*exp(s.nia*e)
 end
diff --git a/src/univariate/continuous/exponential.jl b/src/univariate/continuous/exponential.jl
index d1c487057f..573a469b07 100644
--- a/src/univariate/continuous/exponential.jl
+++ b/src/univariate/continuous/exponential.jl
@@ -105,7 +105,7 @@ 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{T}) where {T} = xval(d, randexp(rng, float(T)))
 
 
 #### Fit model
diff --git a/src/univariate/continuous/tdist.jl b/src/univariate/continuous/tdist.jl
index 4e55e8cfbc..873925c3eb 100644
--- a/src/univariate/continuous/tdist.jl
+++ b/src/univariate/continuous/tdist.jl
@@ -82,7 +82,7 @@ end
 function rand(rng::AbstractRNG, d::TDist)
     ν = d.ν
     z = sqrt(rand(rng, Chisq{typeof(ν)}(ν)) / ν)
-    return randn(rng) / (isinf(ν) ? one(z) : z)
+    return randn(rng, typeof(z)) / (isinf(ν) ? one(z) : z)
 end
 
 function cf(d::TDist{T}, t::Real) where T <: Real
diff --git a/test/univariate/continuous/chisq.jl b/test/univariate/continuous/chisq.jl
index 156bcdb5fa..534b7aaf09 100644
--- a/test/univariate/continuous/chisq.jl
+++ b/test/univariate/continuous/chisq.jl
@@ -1,2 +1,9 @@
-test_cgf(Chisq(1), (0.49, -1, -100, -1f6))
-test_cgf(Chisq(3), (0.49, -1, -100, -1f6))
+
+@testset "Chisq" begin
+    test_cgf(Chisq(1), (0.49, -1, -100, -1.0f6))
+    test_cgf(Chisq(3), (0.49, -1, -100, -1.0f6))
+
+    for T in (Float32, Float64)
+        @test @inferred(rand(Chisq(T(1)))) isa T
+    end
+end
diff --git a/test/univariate/continuous/exponential.jl b/test/univariate/continuous/exponential.jl
index 6704554a79..191528fd7e 100644
--- a/test/univariate/continuous/exponential.jl
+++ b/test/univariate/continuous/exponential.jl
@@ -1,4 +1,10 @@
 
-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))
+@testset "Exponential" begin
+    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))
+
+    for T in (Float32, Float64)
+        @test @inferred(rand(Exponential(T(1)))) isa T
+    end
+end
diff --git a/test/univariate/continuous/gamma.jl b/test/univariate/continuous/gamma.jl
index 3937af191e..8c1887786e 100644
--- a/test/univariate/continuous/gamma.jl
+++ b/test/univariate/continuous/gamma.jl
@@ -1,25 +1,33 @@
 using Test, Distributions, OffsetArrays
 
-test_cgf(Gamma(1  ,1  ), (0.9, -1, -100f0, -1e6))
-test_cgf(Gamma(10 ,1  ), (0.9, -1, -100f0, -1e6))
-test_cgf(Gamma(0.2, 10), (0.08, -1, -100f0, -1e6))
+@testset "Gamma" begin
+    test_cgf(Gamma(1, 1), (0.9, -1, -100.0f0, -1e6))
+    test_cgf(Gamma(10, 1), (0.9, -1, -100.0f0, -1e6))
+    test_cgf(Gamma(0.2, 10), (0.08, -1, -100.0f0, -1e6))
 
-@testset "Gamma suffstats and OffsetArrays" begin
-    a = rand(Gamma(), 11)
-    wa = 1.0:11.0
+    @testset "Gamma suffstats and OffsetArrays" begin
+        a = rand(Gamma(), 11)
+        wa = 1.0:11.0
 
-    resulta = @inferred(suffstats(Gamma, a))
+        resulta = @inferred(suffstats(Gamma, a))
 
-    resultwa = @inferred(suffstats(Gamma, a, wa))
+        resultwa = @inferred(suffstats(Gamma, a, wa))
 
-    b = OffsetArray(a, -5:5)
-    wb = OffsetArray(wa, -5:5)
+        b = OffsetArray(a, -5:5)
+        wb = OffsetArray(wa, -5:5)
 
-    resultb = @inferred(suffstats(Gamma, b))
-    @test resulta == resultb
+        resultb = @inferred(suffstats(Gamma, b))
+        @test resulta == resultb
 
-    resultwb = @inferred(suffstats(Gamma, b, wb))
-    @test resultwa == resultwb
+        resultwb = @inferred(suffstats(Gamma, b, wb))
+        @test resultwa == resultwb
 
-    @test_throws DimensionMismatch suffstats(Gamma, a, wb)
+        @test_throws DimensionMismatch suffstats(Gamma, a, wb)
+    end
+
+    for T in (Float32, Float64)
+        @test @inferred(rand(Gamma(T(1), T(1)))) isa T
+        @test @inferred(rand(Gamma(1/T(2), T(1)))) isa T
+        @test @inferred(rand(Gamma(T(2), T(1)))) isa T
+    end
 end
diff --git a/test/univariate/continuous/tdist.jl b/test/univariate/continuous/tdist.jl
index 16fab2812c..127b992434 100644
--- a/test/univariate/continuous/tdist.jl
+++ b/test/univariate/continuous/tdist.jl
@@ -3,10 +3,17 @@ using ForwardDiff
 
 using Test
 
-@testset "Type stability of `rand` (#1614)" begin
-    if VERSION >= v"1.9.0-DEV.348"
-        # randn(::BigFloat) was only added in https://github.com/JuliaLang/julia/pull/44714
-        @inferred(rand(TDist(big"1.0")))
+@testset "TDist" begin
+    @testset "Type stability of `rand` (#1614)" begin
+        if VERSION >= v"1.9.0-DEV.348"
+            # randn(::BigFloat) was only added in https://github.com/JuliaLang/julia/pull/44714
+            @inferred(rand(TDist(big"1.0")))
+        end
+        @inferred(rand(TDist(ForwardDiff.Dual(1.0))))
+
+    end
+
+    for T in (Float32, Float64)
+        @test @inferred(rand(TDist(T(1)))) isa T
     end
-    @inferred(rand(TDist(ForwardDiff.Dual(1.0))))
 end

From 00b7fad421c606c1f7d58c38ab0114472e76a747 Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Fri, 23 Aug 2024 21:44:39 +0200
Subject: [PATCH 13/19] Update Project.toml

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

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

From aad76a3a8d559bd89391d67d9de98874f89d9744 Mon Sep 17 00:00:00 2001
From: David Widmann <devmotion@users.noreply.github.com>
Date: Mon, 26 Aug 2024 17:02:24 +0200
Subject: [PATCH 14/19] Test pre-release Julia versions (#1888)

---
 .github/workflows/CI.yml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 384a02a315..886724bd29 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -25,6 +25,7 @@ jobs:
         version:
           - '1.3'
           - '1'
+          - pre
         os:
           - ubuntu-latest
           - macos-latest

From e1340f02b5604cedf02df93b9ae202f3679d9818 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Mon, 26 Aug 2024 23:51:43 +0530
Subject: [PATCH 15/19] Relax Fill return type in mvlognormal tests (#1890)

* Relax Fill return type in mvlognormal tests

* Use AbstractFill instead of Union

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

---------

Co-authored-by: David Widmann <devmotion@users.noreply.github.com>
---
 test/multivariate/mvlognormal.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/multivariate/mvlognormal.jl b/test/multivariate/mvlognormal.jl
index 89f923b889..7ef76bef2d 100644
--- a/test/multivariate/mvlognormal.jl
+++ b/test/multivariate/mvlognormal.jl
@@ -19,7 +19,7 @@ function test_mvlognormal(g::MvLogNormal, n_tsamples::Int=10^6,
     @test partype(g) == Float64
     @test isa(mn, Vector{Float64})
     if g.normal.μ isa Zeros{Float64,1}
-        @test md isa Fill{Float64,1}
+        @test md isa FillArrays.AbstractFill{Float64,1}
     else
         @test md isa Vector{Float64}
     end

From b219803a0d03a7c75d7aef7c0bab6cd0d79997dc Mon Sep 17 00:00:00 2001
From: "Steven G. Johnson" <stevenj@mit.edu>
Date: Mon, 2 Sep 2024 02:45:02 -0400
Subject: [PATCH 16/19] mark the MatrixReshaped type as deprecated (#1893)

* mark the MatrixReshaped type as deprecated

* Update deprecates.jl
---
 src/deprecates.jl | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/deprecates.jl b/src/deprecates.jl
index 3cfc20a5c6..3106154a5e 100644
--- a/src/deprecates.jl
+++ b/src/deprecates.jl
@@ -54,6 +54,7 @@ end
 
 # Deprecate `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))

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 17/19] 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 18/19] 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 19/19] 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"