Merge branch 'master' into namedtuplevariate

.github/workflows/CI.yml

@@ -25,6 +25,7 @@ jobs:
         version:
           - '1.3'
           - '1'
+          - pre
         os:
           - ubuntu-latest
           - macos-latest

Project.toml

@@ -1,7 +1,7 @@
 name = "Distributions"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
 authors = ["JuliaStats"]
-version = "0.25.109"
+version = "0.25.111"
 
 [deps]
 AliasTables = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"

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))

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

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

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

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

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

test/runtests.jl

@@ -97,6 +97,7 @@ const tests = [
     "eachvariate",
     "univariate/continuous/triangular",
     "statsapi",
+    "univariate/continuous/inversegaussian",
 
     ### missing files compared to /src:
     # "common",
@@ -138,7 +139,6 @@ const tests = [
     # "univariate/continuous/generalizedextremevalue",
     # "univariate/continuous/generalizedpareto",
     # "univariate/continuous/inversegamma",
-    # "univariate/continuous/inversegaussian",
     # "univariate/continuous/ksdist",
     # "univariate/continuous/ksonesided",
     # "univariate/continuous/levy",

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

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

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

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

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