Merge branch 'master' into roots.jl

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"

src/Distributions.jl

@@ -129,7 +129,6 @@ export
     MatrixBeta,
     MatrixFDist,
     MatrixNormal,
-    MatrixReshaped,
     MatrixTDist,
     MixtureModel,
     Multinomial,

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)

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
 

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
 

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
 

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
 

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
 

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
 

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

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

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

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