From ef35e00137c1381e25cdc3516940e54e57b0746b Mon Sep 17 00:00:00 2001 From: David Widmann Date: Sat, 22 Jan 2022 01:41:17 +0100 Subject: [PATCH 1/7] Use Julia implementations of `normpdf` etc. in StatsFuns --- Project.toml | 2 +- src/univariate/continuous/normal.jl | 159 +--------------------------- 2 files changed, 4 insertions(+), 157 deletions(-) diff --git a/Project.toml b/Project.toml index 598facf61..82f507ccc 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Distributions" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" authors = ["JuliaStats"] -version = "0.25.41" +version = "0.25.42" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" diff --git a/src/univariate/continuous/normal.jl b/src/univariate/continuous/normal.jl index db896ced8..89f461e3f 100644 --- a/src/univariate/continuous/normal.jl +++ b/src/univariate/continuous/normal.jl @@ -85,163 +85,10 @@ end #### Evaluation -# Helpers -""" - xval(d::Normal, z::Real) - -Computes the x-value based on a Normal distribution and a z-value. -""" -function xval(d::Normal, z::Real) - if isinf(z) && iszero(d.σ) - d.μ + one(d.σ) * z - else - d.μ + d.σ * z - end -end -""" - zval(d::Normal, x::Real) - -Computes the z-value based on a Normal distribution and a x-value. -""" -zval(d::Normal, x::Real) = (x - d.μ) / d.σ - -gradlogpdf(d::Normal, x::Real) = -zval(d, x) / d.σ - -# logpdf -_normlogpdf(z::Real) = -(abs2(z) + log2π)/2 - -function logpdf(d::Normal, x::Real) - μ, σ = d.μ, d.σ - if iszero(d.σ) - if x == μ - z = zval(Normal(μ, one(σ)), x) - else - z = zval(d, x) - σ = one(σ) - end - else - z = zval(Normal(μ, σ), x) - end - return _normlogpdf(z) - log(σ) -end +# Use Julia implementations in StatsFuns +@_delegate_statsfuns Normal norm μ σ -# pdf -_normpdf(z::Real) = exp(-abs2(z)/2) * invsqrt2π - -function pdf(d::Normal, x::Real) - μ, σ = d.μ, d.σ - if iszero(σ) - if x == μ - z = zval(Normal(μ, one(σ)), x) - else - z = zval(d, x) - σ = one(σ) - end - else - z = zval(Normal(μ, σ), x) - end - return _normpdf(z) / σ -end - -# logcdf -function _normlogcdf(z::Real) - if z < -one(z) - return log(erfcx(-z * invsqrt2)/2) - abs2(z)/2 - else - return log1p(-erfc(z * invsqrt2)/2) - end -end - -function logcdf(d::Normal, x::Real) - if iszero(d.σ) && x == d.μ - z = zval(Normal(zero(d.μ), d.σ), one(x)) - else - z = zval(d, x) - end - return _normlogcdf(z) -end - -# logccdf -function _normlogccdf(z::Real) - if z > one(z) - return log(erfcx(z * invsqrt2)/2) - abs2(z)/2 - else - return log1p(-erfc(-z * invsqrt2)/2) - end -end - -function logccdf(d::Normal, x::Real) - if iszero(d.σ) && x == d.μ - z = zval(Normal(zero(d.μ), d.σ), one(x)) - else - z = zval(d, x) - end - return _normlogccdf(z) -end - -# cdf -_normcdf(z::Real) = erfc(-z * invsqrt2)/2 - -function cdf(d::Normal, x::Real) - if iszero(d.σ) && x == d.μ - z = zval(Normal(zero(d.μ), d.σ), one(x)) - else - z = zval(d, x) - end - return _normcdf(z) -end - -# ccdf -_normccdf(z::Real) = erfc(z * invsqrt2)/2 - -function ccdf(d::Normal, x::Real) - if iszero(d.σ) && x == d.μ - z = zval(Normal(zero(d.μ), d.σ), one(x)) - else - z = zval(d, x) - end - return _normccdf(z) -end - -# quantile -function quantile(d::Normal, p::Real) - # Promote to ensure that we don't compute erfcinv in low precision and then promote - _p, _μ, _σ = map(float, promote(p, d.μ, d.σ)) - q = xval(d, -erfcinv(2*_p) * sqrt2) - if isnan(_p) - return oftype(q, _p) - elseif iszero(_σ) - # Quantile not uniquely defined at p=0 and p=1 when σ=0 - if iszero(_p) - return oftype(q, -Inf) - elseif isone(_p) - return oftype(q, Inf) - else - return oftype(q, _μ) - end - end - return q -end - -# cquantile -function cquantile(d::Normal, p::Real) - # Promote to ensure that we don't compute erfcinv in low precision and then promote - _p, _μ, _σ = map(float, promote(p, d.μ, d.σ)) - q = xval(d, erfcinv(2*_p) * sqrt2) - if isnan(_p) - return oftype(q, _p) - elseif iszero(_σ) - # Quantile not uniquely defined at p=0 and p=1 when σ=0 - if iszero(_p) - return oftype(q, Inf) - elseif isone(_p) - return oftype(q, -Inf) - else - return oftype(q, _μ) - end - end - return q -end +gradlogpdf(d::Normal, x::Real) = -(x - d.μ) / d.σ^2 mgf(d::Normal, t::Real) = exp(t * d.μ + d.σ^2 / 2 * t^2) cf(d::Normal, t::Real) = exp(im * t * d.μ - d.σ^2 / 2 * t^2) From 7833c62b042e35b15bbe12026dad27197698f262 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Mon, 24 Jan 2022 23:20:14 +0100 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 82f507ccc..75d494a7f 100644 --- a/Project.toml +++ b/Project.toml @@ -27,7 +27,7 @@ PDMats = "0.10, 0.11" QuadGK = "2" SpecialFunctions = "0.8, 0.9, 0.10, 1.0, 2" StatsBase = "0.32, 0.33" -StatsFuns = "0.8, 0.9" +StatsFuns = "0.9.15" julia = "1" [extras] From bc2c60dced6e2babacb45176ed1aab3f6ecc8784 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Tue, 25 Jan 2022 15:07:42 +0100 Subject: [PATCH 3/7] Skip tests on Julia < 1.3 --- test/normal.jl | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/test/normal.jl b/test/normal.jl index 487dcebba..7c3bc3b0f 100644 --- a/test/normal.jl +++ b/test/normal.jl @@ -149,7 +149,11 @@ end @test @inferred(quantile(Normal(1.0f0, 0.0f0), 1.0f0)) === Inf32 @test @inferred(quantile(Normal(1.0f0, 0.0f0), 0.5f0)) === 1.0f0 @test isnan_type(Float32, @inferred(quantile(Normal(1.0f0, 0.0f0), NaN32))) - @test @inferred(quantile(Normal(1//1, 0//1), 1//2)) === 1.0 + # `erfcinv(::Rational)` is defined only in SpecialFunctions >= 1.2.0 which requires Julia >= 1.3 + # Ref: https://github.com/JuliaStats/Distributions.jl/pull/1487#issuecomment-1020626771 + if VERSION >= v"1.3" + @test @inferred(quantile(Normal(1//1, 0//1), 1//2)) === 1.0 + end @test @inferred(quantile(Normal(1f0, 0f0), 1//2)) === 1f0 @test @inferred(quantile(Normal(1f0, 0.0), 1//2)) === 1.0 @@ -161,7 +165,11 @@ end @test @inferred(cquantile(Normal(1.0f0, 0.0f0), 1.0f0)) === -Inf32 @test @inferred(cquantile(Normal(1.0f0, 0.0f0), 0.5f0)) === 1.0f0 @test isnan_type(Float32, @inferred(cquantile(Normal(1.0f0, 0.0f0), NaN32))) - @test @inferred(cquantile(Normal(1//1, 0//1), 1//2)) === 1.0 + # `erfcinv(::Rational)` is defined only in SpecialFunctions >= 1.2.0 which requires Julia >= 1.3 + # Ref: https://github.com/JuliaStats/Distributions.jl/pull/1487#issuecomment-1020626771 + if VERSION >= v"1.3" + @test @inferred(cquantile(Normal(1//1, 0//1), 1//2)) === 1.0 + end @test @inferred(cquantile(Normal(1f0, 0f0), 1//2)) === 1f0 @test @inferred(cquantile(Normal(1f0, 0.0), 1//2)) === 1.0 end From 0c19af21058591ef1cf6a6043fd83fb34f0d20d0 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Tue, 25 Jan 2022 15:11:04 +0100 Subject: [PATCH 4/7] Simplify `gradlogpdf` --- src/univariate/continuous/normal.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/univariate/continuous/normal.jl b/src/univariate/continuous/normal.jl index 89f461e3f..5b42dd5c2 100644 --- a/src/univariate/continuous/normal.jl +++ b/src/univariate/continuous/normal.jl @@ -88,7 +88,7 @@ end # Use Julia implementations in StatsFuns @_delegate_statsfuns Normal norm μ σ -gradlogpdf(d::Normal, x::Real) = -(x - d.μ) / d.σ^2 +gradlogpdf(d::Normal, x::Real) = (d.μ - x) / d.σ^2 mgf(d::Normal, t::Real) = exp(t * d.μ + d.σ^2 / 2 * t^2) cf(d::Normal, t::Real) = exp(im * t * d.μ - d.σ^2 / 2 * t^2) From 529b9071e2280dc21d3de2a31755fca4062d4df8 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Tue, 25 Jan 2022 15:28:15 +0100 Subject: [PATCH 5/7] Fix `LogNormal` tests --- test/lognormal.jl | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/test/lognormal.jl b/test/lognormal.jl index 0806356b7..4532555a0 100644 --- a/test/lognormal.jl +++ b/test/lognormal.jl @@ -265,7 +265,11 @@ end @test @inferred(quantile(LogNormal(1.0f0, 0.0f0), 1.0f0)) === Inf32 @test @inferred(quantile(LogNormal(1.0f0, 0.0f0), 0.5f0)) === exp(1.0f0) @test isnan_type(Float32, @inferred(quantile(LogNormal(1.0f0, 0.0f0), NaN32))) - @test @inferred(quantile(LogNormal(1//1, 0//1), 1//2)) === exp(1) + # `erfcinv(::Rational)` is defined only in SpecialFunctions >= 1.2.0 which requires Julia >= 1.3 + # Ref: https://github.com/JuliaStats/Distributions.jl/pull/1487#issuecomment-1020626771 + if VERSION >= v"1.3" + @test @inferred(quantile(LogNormal(1//1, 0//1), 1//2)) === exp(1) + end # cquantile @test @inferred(cquantile(LogNormal(1.0, 0.0), 0.0f0)) === Inf @@ -276,7 +280,11 @@ end @test @inferred(cquantile(LogNormal(1.0f0, 0.0f0), 1.0f0)) === 0.0f0 @test @inferred(cquantile(LogNormal(1.0f0, 0.0f0), 0.5f0)) === exp(1.0f0) @test isnan_type(Float32, @inferred(cquantile(LogNormal(1.0f0, 0.0f0), NaN32))) - @test @inferred(cquantile(LogNormal(1//1, 0//1), 1//2)) === exp(1) + # `erfcinv(::Rational)` is defined only in SpecialFunctions >= 1.2.0 which requires Julia >= 1.3 + # Ref: https://github.com/JuliaStats/Distributions.jl/pull/1487#issuecomment-1020626771 + if VERSION >= v"1.3" + @test @inferred(cquantile(LogNormal(1//1, 0//1), 1//2)) === exp(1) + end # gradlogpdf @test @inferred(gradlogpdf(LogNormal(0.0, 1.0), 1.0)) === -1.0 From cd880c8edf8bef92a236f9a18494d8a4768baf60 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Tue, 1 Feb 2022 00:59:45 +0100 Subject: [PATCH 6/7] Require Julia >= 1.3 and SpecialFunctions >= 1.2 --- .github/workflows/CI.yml | 2 +- Project.toml | 4 ++-- test/lognormal.jl | 16 ++++------------ test/normal.jl | 16 ++++------------ 4 files changed, 11 insertions(+), 27 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 9d478dc99..cfa5a4762 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -24,7 +24,7 @@ jobs: fail-fast: false matrix: version: - - '1.0' + - '1.3' - '1' - 'nightly' os: diff --git a/Project.toml b/Project.toml index 75d494a7f..8eee56e1f 100644 --- a/Project.toml +++ b/Project.toml @@ -25,10 +25,10 @@ DensityInterface = "0.4" FillArrays = "0.9, 0.10, 0.11, 0.12" PDMats = "0.10, 0.11" QuadGK = "2" -SpecialFunctions = "0.8, 0.9, 0.10, 1.0, 2" +SpecialFunctions = "1.2, 2" StatsBase = "0.32, 0.33" StatsFuns = "0.9.15" -julia = "1" +julia = "1.3" [extras] Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" diff --git a/test/lognormal.jl b/test/lognormal.jl index 4532555a0..f43b99188 100644 --- a/test/lognormal.jl +++ b/test/lognormal.jl @@ -17,8 +17,8 @@ isnan_type(::Type{T}, v) where {T} = isnan(v) && v isa T @test logdiffcdf(LogNormal(), Float64(exp(5)), Float64(exp(3))) ≈ -6.607938594596893 rtol=1e-12 let d = LogNormal(Float64(0), Float64(1)), x = Float64(exp(-60)), y = Float64(exp(-60.001)) float_res = logdiffcdf(d, x, y) - big_x = VERSION < v"1.1" ? BigFloat(x, 100) : BigFloat(x; precision=100) - big_y = VERSION < v"1.1" ? BigFloat(y, 100) : BigFloat(y; precision=100) + big_x = BigFloat(x; precision=100) + big_y = BigFloat(y; precision=100) big_float_res = log(cdf(d, big_x) - cdf(d, big_y)) @test float_res ≈ big_float_res end @@ -265,11 +265,7 @@ end @test @inferred(quantile(LogNormal(1.0f0, 0.0f0), 1.0f0)) === Inf32 @test @inferred(quantile(LogNormal(1.0f0, 0.0f0), 0.5f0)) === exp(1.0f0) @test isnan_type(Float32, @inferred(quantile(LogNormal(1.0f0, 0.0f0), NaN32))) - # `erfcinv(::Rational)` is defined only in SpecialFunctions >= 1.2.0 which requires Julia >= 1.3 - # Ref: https://github.com/JuliaStats/Distributions.jl/pull/1487#issuecomment-1020626771 - if VERSION >= v"1.3" - @test @inferred(quantile(LogNormal(1//1, 0//1), 1//2)) === exp(1) - end + @test @inferred(quantile(LogNormal(1//1, 0//1), 1//2)) === exp(1) # cquantile @test @inferred(cquantile(LogNormal(1.0, 0.0), 0.0f0)) === Inf @@ -280,11 +276,7 @@ end @test @inferred(cquantile(LogNormal(1.0f0, 0.0f0), 1.0f0)) === 0.0f0 @test @inferred(cquantile(LogNormal(1.0f0, 0.0f0), 0.5f0)) === exp(1.0f0) @test isnan_type(Float32, @inferred(cquantile(LogNormal(1.0f0, 0.0f0), NaN32))) - # `erfcinv(::Rational)` is defined only in SpecialFunctions >= 1.2.0 which requires Julia >= 1.3 - # Ref: https://github.com/JuliaStats/Distributions.jl/pull/1487#issuecomment-1020626771 - if VERSION >= v"1.3" - @test @inferred(cquantile(LogNormal(1//1, 0//1), 1//2)) === exp(1) - end + @test @inferred(cquantile(LogNormal(1//1, 0//1), 1//2)) === exp(1) # gradlogpdf @test @inferred(gradlogpdf(LogNormal(0.0, 1.0), 1.0)) === -1.0 diff --git a/test/normal.jl b/test/normal.jl index 7c3bc3b0f..7682c17e7 100644 --- a/test/normal.jl +++ b/test/normal.jl @@ -14,8 +14,8 @@ isnan_type(::Type{T}, v) where {T} = isnan(v) && v isa T @test logdiffcdf(Normal(), Float64(5), Float64(3)) ≈ -6.607938594596893 rtol=1e-12 let d = Normal(Float64(0), Float64(1)), x = Float64(-60), y = Float64(-60.001) float_res = logdiffcdf(d, x, y) - big_x = VERSION < v"1.1" ? BigFloat(x, 100) : BigFloat(x; precision=100) - big_y = VERSION < v"1.1" ? BigFloat(y, 100) : BigFloat(y; precision=100) + big_x = BigFloat(x; precision=100) + big_y = BigFloat(y; precision=100) big_float_res = log(cdf(d, big_x) - cdf(d, big_y)) @test float_res ≈ big_float_res end @@ -149,11 +149,7 @@ end @test @inferred(quantile(Normal(1.0f0, 0.0f0), 1.0f0)) === Inf32 @test @inferred(quantile(Normal(1.0f0, 0.0f0), 0.5f0)) === 1.0f0 @test isnan_type(Float32, @inferred(quantile(Normal(1.0f0, 0.0f0), NaN32))) - # `erfcinv(::Rational)` is defined only in SpecialFunctions >= 1.2.0 which requires Julia >= 1.3 - # Ref: https://github.com/JuliaStats/Distributions.jl/pull/1487#issuecomment-1020626771 - if VERSION >= v"1.3" - @test @inferred(quantile(Normal(1//1, 0//1), 1//2)) === 1.0 - end + @test @inferred(quantile(Normal(1//1, 0//1), 1//2)) === 1.0 @test @inferred(quantile(Normal(1f0, 0f0), 1//2)) === 1f0 @test @inferred(quantile(Normal(1f0, 0.0), 1//2)) === 1.0 @@ -165,11 +161,7 @@ end @test @inferred(cquantile(Normal(1.0f0, 0.0f0), 1.0f0)) === -Inf32 @test @inferred(cquantile(Normal(1.0f0, 0.0f0), 0.5f0)) === 1.0f0 @test isnan_type(Float32, @inferred(cquantile(Normal(1.0f0, 0.0f0), NaN32))) - # `erfcinv(::Rational)` is defined only in SpecialFunctions >= 1.2.0 which requires Julia >= 1.3 - # Ref: https://github.com/JuliaStats/Distributions.jl/pull/1487#issuecomment-1020626771 - if VERSION >= v"1.3" - @test @inferred(cquantile(Normal(1//1, 0//1), 1//2)) === 1.0 - end + @test @inferred(cquantile(Normal(1//1, 0//1), 1//2)) === 1.0 @test @inferred(cquantile(Normal(1f0, 0f0), 1//2)) === 1f0 @test @inferred(cquantile(Normal(1f0, 0.0), 1//2)) === 1.0 end From 86113fd840ba81d9cd26ad3b5b65b4911db13d16 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Tue, 1 Feb 2022 01:03:29 +0100 Subject: [PATCH 7/7] Revert changes --- test/normal.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/normal.jl b/test/normal.jl index 7682c17e7..28dbd7a2c 100644 --- a/test/normal.jl +++ b/test/normal.jl @@ -149,7 +149,7 @@ end @test @inferred(quantile(Normal(1.0f0, 0.0f0), 1.0f0)) === Inf32 @test @inferred(quantile(Normal(1.0f0, 0.0f0), 0.5f0)) === 1.0f0 @test isnan_type(Float32, @inferred(quantile(Normal(1.0f0, 0.0f0), NaN32))) - @test @inferred(quantile(Normal(1//1, 0//1), 1//2)) === 1.0 + @test @inferred(quantile(Normal(1//1, 0//1), 1//2)) === 1.0 @test @inferred(quantile(Normal(1f0, 0f0), 1//2)) === 1f0 @test @inferred(quantile(Normal(1f0, 0.0), 1//2)) === 1.0 @@ -161,7 +161,7 @@ end @test @inferred(cquantile(Normal(1.0f0, 0.0f0), 1.0f0)) === -Inf32 @test @inferred(cquantile(Normal(1.0f0, 0.0f0), 0.5f0)) === 1.0f0 @test isnan_type(Float32, @inferred(cquantile(Normal(1.0f0, 0.0f0), NaN32))) - @test @inferred(cquantile(Normal(1//1, 0//1), 1//2)) === 1.0 + @test @inferred(cquantile(Normal(1//1, 0//1), 1//2)) === 1.0 @test @inferred(cquantile(Normal(1f0, 0f0), 1//2)) === 1f0 @test @inferred(cquantile(Normal(1f0, 0.0), 1//2)) === 1.0 end