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 1/7] Fix InverseGaussian cdf overflows (#1882) * Fix inversegaussian cdf overflows * Update src/univariate/continuous/inversegaussian.jl Co-authored-by: David Widmann * Update src/univariate/continuous/inversegaussian.jl Co-authored-by: David Widmann * Fix missing parentheses * Move parentheses * Move new tests to dedicated inversegaussian file * Add ccdf tests --------- Co-authored-by: David Widmann --- 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 Date: Fri, 9 Aug 2024 11:38:22 +0200 Subject: [PATCH 2/7] 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 Date: Fri, 23 Aug 2024 15:44:22 -0400 Subject: [PATCH 3/7] 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 * fix type specificaton in `rand(TDist)` Co-authored-by: David Widmann * 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 * fix use `@inferred` in tests for `Gamma` Co-authored-by: David Widmann * fix use `@inferred` in tests for `TDist` Co-authored-by: David Widmann * 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 --- 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 Date: Fri, 23 Aug 2024 21:44:39 +0200 Subject: [PATCH 4/7] 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 Date: Mon, 26 Aug 2024 17:02:24 +0200 Subject: [PATCH 5/7] 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 Date: Mon, 26 Aug 2024 23:51:43 +0530 Subject: [PATCH 6/7] 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 --------- Co-authored-by: David Widmann --- 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" Date: Mon, 2 Sep 2024 02:45:02 -0400 Subject: [PATCH 7/7] 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))