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 1e8a71b70..b581f9d0b 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Distributions" uuid = "31c24e10-a181-5473-b8eb-7969acd0382f" authors = ["JuliaStats"] -version = "0.25.45" +version = "0.25.46" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" @@ -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.8, 0.9" -julia = "1" +StatsFuns = "0.9.15" +julia = "1.3" [extras] Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9" diff --git a/src/univariate/continuous/normal.jl b/src/univariate/continuous/normal.jl index 93043fd60..d21e95d8e 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) = (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) diff --git a/test/lognormal.jl b/test/lognormal.jl index 0806356b7..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 diff --git a/test/normal.jl b/test/normal.jl index 487dcebba..28dbd7a2c 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