Use Julia implementations of normpdf etc. in StatsFuns

.github/workflows/CI.yml

@@ -24,7 +24,7 @@ jobs:
  fail-fast: false
  matrix:
  version:
- - '1.0'
+ - '1.3'
  - '1'
  - 'nightly'
  os:

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"

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)

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

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