Skip to content

Commit 5eda5ca

Browse files
devmotiondevmotion
committed
Revert "Revert "Use julia implementations for pdfs and some cdf-like functions (#113)""
This reverts commit 6220584.
1 parent cede7a9 commit 5eda5ca

File tree

9 files changed

+278
-95
lines changed

9 files changed

+278
-95
lines changed

‎Project.toml‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@ version = "0.9.17"
44

55
[deps]
66
ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
7+
HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
78
InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112"
89
IrrationalConstants = "92d709cd-6900-40b7-9082-c6be49f344b6"
910
LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
@@ -13,6 +14,7 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
1314

1415
[compat]
1516
ChainRulesCore = "1"
17+
HypergeometricFunctions = "0.3"
1618
InverseFunctions = "0.1"
1719
IrrationalConstants = "0.1"
1820
LogExpFunctions = "0.3.2"

‎src/distrs/beta.jl‎

Lines changed: 67 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,17 @@
11
# functions related to beta distributions
22

3+
using HypergeometricFunctions: _₂F₁
4+
35
# R implementations
4-
# For pdf and logpdf we use the Julia implementation
56
using .RFunctions:
6-
betacdf,
7-
betaccdf,
8-
betalogcdf,
9-
betalogccdf,
10-
betainvcdf,
11-
betainvccdf,
7+
# betapdf,
8+
# betalogpdf,
9+
# betacdf,
10+
# betaccdf,
11+
# betalogcdf,
12+
# betalogccdf,
13+
# betainvcdf,
14+
# betainvccdf,
1215
betainvlogcdf,
1316
betainvlogccdf
1417

@@ -22,3 +25,60 @@ function betalogpdf(α::T, β::T, x::T) where {T<:Real}
2225
val = xlogy(α - 1, y) + xlog1py(β - 1, -y) - logbeta(α, β)
2326
return x < 0 || x > 1 ? oftype(val, -Inf) : val
2427
end
28+
29+
function betacdf(α::Real, β::Real, x::Real)
30+
# Handle a degenerate case
31+
if iszero(α) && β > 0
32+
return float(last(promote(α, β, x, x >= 0)))
33+
end
34+
35+
return first(beta_inc(α, β, clamp(x, 0, 1)))
36+
end
37+
38+
function betaccdf(α::Real, β::Real, x::Real)
39+
# Handle a degenerate case
40+
if iszero(α) && β > 0
41+
return float(last(promote(α, β, x, x < 0)))
42+
end
43+
44+
last(beta_inc(α, β, clamp(x, 0, 1)))
45+
end
46+
47+
# The log version is currently based on non-log version. When the cdf is very small we shift
48+
# to an implementation based on the hypergeometric function ₂F₁ to avoid underflow.
49+
function betalogcdf(α::T, β::T, x::T) where {T<:Real}
50+
# Handle a degenerate case
51+
if iszero(α) && β > 0
52+
return log(last(promote(x, x >= 0)))
53+
end
54+
55+
_x = clamp(x, 0, 1)
56+
p, q = beta_inc(α, β, _x)
57+
if p < floatmin(p)
58+
# see https://dlmf.nist.gov/8.17#E7
59+
return -log(α) + xlogy(α, _x) + log(_₂F₁(promote(α, 1 - β, α + 1, _x)...)) - logbeta(α, β)
60+
elseif p <= 0.7
61+
return log(p)
62+
else
63+
return log1p(-q)
64+
end
65+
end
66+
betalogcdf(α::Real, β::Real, x::Real) = betalogcdf(promote(α, β, x)...)
67+
68+
function betalogccdf(α::Real, β::Real, x::Real)
69+
# Handle a degenerate case
70+
if iszero(α) && β > 0
71+
return log(last(promote(α, β, x, x < 0)))
72+
end
73+
74+
p, q = beta_inc(α, β, clamp(x, 0, 1))
75+
if q < 0.7
76+
return log(q)
77+
else
78+
return log1p(-p)
79+
end
80+
end
81+
82+
betainvcdf(α::Real, β::Real, p::Real) = first(beta_inc_inv(α, β, p))
83+
84+
betainvccdf(α::Real, β::Real, p::Real) = last(beta_inc_inv(β, α, p))

‎src/distrs/binom.jl‎

Lines changed: 26 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,18 +1,18 @@
11
# functions related to binomial distribution
22

33
# R implementations
4-
# For pdf and logpdf we use the Julia implementation
54
using .RFunctions:
6-
binomcdf,
7-
binomccdf,
8-
binomlogcdf,
9-
binomlogccdf,
5+
# binompdf,
6+
# binomlogpdf,
7+
# binomcdf,
8+
# binomccdf,
9+
# binomlogcdf,
10+
# binomlogccdf,
1011
binominvcdf,
1112
binominvccdf,
1213
binominvlogcdf,
1314
binominvlogccdf
1415

15-
1616
# Julia implementations
1717
binompdf(n::Real, p::Real, k::Real) = exp(binomlogpdf(n, p, k))
1818

@@ -22,3 +22,23 @@ function binomlogpdf(n::T, p::T, k::T) where {T<:Real}
2222
val = min(0, betalogpdf(m + 1, n - m + 1, p) - log(n + 1))
2323
return 0 <= k <= n && isinteger(k) ? val : oftype(val, -Inf)
2424
end
25+
26+
for l in ("", "log"), compl in (false, true)
27+
fbinom = Symbol(string("binom", l, ifelse(compl, "c", "" ), "cdf"))
28+
fbeta = Symbol(string("beta" , l, ifelse(compl, "", "c"), "cdf"))
29+
@eval function ($fbinom)(n::Real, p::Real, k::Real)
30+
if isnan(k)
31+
return last(promote(n, p, k))
32+
end
33+
res = ($fbeta)(max(0, floor(k) + 1), max(0, n - floor(k)), p)
34+
35+
# When p == 1 the betaccdf doesn't return the correct result
36+
# so these cases have to be special cased
37+
if isone(p)
38+
newres = oftype(res, $compl ? k < n : k >= n)
39+
return $(l === "" ? :newres : :(log(newres)))
40+
else
41+
return res
42+
end
43+
end
44+
end

‎src/distrs/chisq.jl‎

Lines changed: 9 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -1,22 +1,11 @@
11
# functions related to chi-square distribution
22

3-
# R implementations
4-
# For pdf and logpdf we use the Julia implementation
5-
using .RFunctions:
6-
chisqcdf,
7-
chisqccdf,
8-
chisqlogcdf,
9-
chisqlogccdf,
10-
chisqinvcdf,
11-
chisqinvccdf,
12-
chisqinvlogcdf,
13-
chisqinvlogccdf
14-
15-
# Julia implementations
16-
# promotion ensures that we do forward e.g. `chisqpdf(::Int, ::Float32)` to
17-
# `gammapdf(::Float32, ::Int, ::Float32)` but not `gammapdf(::Float64, ::Int, ::Float32)`
18-
chisqpdf(k::Real, x::Real) = chisqpdf(promote(k, x)...)
19-
chisqpdf(k::T, x::T) where {T<:Real} = gammapdf(k / 2, 2, x)
20-
21-
chisqlogpdf(k::Real, x::Real) = chisqlogpdf(promote(k, x)...)
22-
chisqlogpdf(k::T, x::T) where {T<:Real} = gammalogpdf(k / 2, 2, x)
3+
# Just use the Gamma definitions
4+
for f in ("pdf", "logpdf", "cdf", "ccdf", "logcdf", "logccdf", "invcdf", "invccdf", "invlogcdf", "invlogccdf")
5+
_chisqf = Symbol("chisq"*f)
6+
_gammaf = Symbol("gamma"*f)
7+
@eval begin
8+
$(_chisqf)(k::Real, x::Real) = $(_chisqf)(promote(k, x)...)
9+
$(_chisqf)(k::T, x::T) where {T<:Real} = $(_gammaf)(k/2, 2, x)
10+
end
11+
end

‎src/distrs/fdist.jl‎

Lines changed: 16 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,5 @@
11
# functions related to F distribution
22

3-
# R implementations
4-
# For pdf and logpdf we use the Julia implementation
5-
using .RFunctions:
6-
fdistcdf,
7-
fdistccdf,
8-
fdistlogcdf,
9-
fdistlogccdf,
10-
fdistinvcdf,
11-
fdistinvccdf,
12-
fdistinvlogcdf,
13-
fdistinvlogccdf
14-
153
# Julia implementations
164
fdistpdf(ν1::Real, ν2::Real, x::Real) = exp(fdistlogpdf(ν1, ν2, x))
175

@@ -23,3 +11,19 @@ function fdistlogpdf(ν1::T, ν2::T, x::T) where {T<:Real}
2311
val = (xlogy(ν1, ν1ν2) + xlogy(ν1 - 2, y) - xlogy(ν1 + ν2, 1 + ν1ν2 * y)) / 2 - logbeta(ν1 / 2, ν2 / 2)
2412
return x < 0 ? oftype(val, -Inf) : val
2513
end
14+
15+
for f in ("cdf", "ccdf", "logcdf", "logccdf")
16+
ff = Symbol("fdist"*f)
17+
bf = Symbol("beta"*f)
18+
@eval $ff(ν1::T, ν2::T, x::T) where {T<:Real} = $bf(ν1/2, ν2/2, inv(1 + ν2/(ν1*max(0, x))))
19+
@eval $ff(ν1::Real, ν2::Real, x::Real) = $ff(promote(ν1, ν2, x)...)
20+
end
21+
for f in ("invcdf", "invccdf", "invlogcdf", "invlogccdf")
22+
ff = Symbol("fdist"*f)
23+
bf = Symbol("beta"*f)
24+
@eval function $ff(ν1::T, ν2::T, y::T) where {T<:Real}
25+
x = $bf(ν1/2, ν2/2, y)
26+
return x/(1 - x)*ν2/ν1
27+
end
28+
@eval $ff(ν1::Real, ν2::Real, y::Real) = $ff(promote(ν1, ν2, y)...)
29+
end

‎src/distrs/gamma.jl‎

Lines changed: 87 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,17 @@
11
# functions related to gamma distribution
22

3+
using HypergeometricFunctions: drummond1F1
4+
35
# R implementations
4-
# For pdf and logpdf we use the Julia implementation
56
using .RFunctions:
6-
gammacdf,
7-
gammaccdf,
8-
gammalogcdf,
9-
gammalogccdf,
10-
gammainvcdf,
11-
gammainvccdf,
7+
# gammapdf,
8+
# gammalogpdf,
9+
# gammacdf,
10+
# gammaccdf,
11+
# gammalogcdf,
12+
# gammalogccdf,
13+
# gammainvcdf,
14+
# gammainvccdf,
1215
gammainvlogcdf,
1316
gammainvlogccdf
1417

@@ -22,3 +25,80 @@ function gammalogpdf(k::T, θ::T, x::T) where {T<:Real}
2225
val = -loggamma(k) + xlogy(k - 1, xθ) - log(θ) - xθ
2326
return x < 0 ? oftype(val, -Inf) : val
2427
end
28+
29+
function gammacdf(k::T, θ::T, x::T) where {T<:Real}
30+
# Handle the degenerate case
31+
if iszero(k)
32+
return float(last(promote(x, x >= 0)))
33+
end
34+
return first(gamma_inc(k, max(0, x)/θ))
35+
end
36+
gammacdf(k::Real, θ::Real, x::Real) = gammacdf(map(float, promote(k, θ, x))...)
37+
38+
function gammaccdf(k::T, θ::T, x::T) where {T<:Real}
39+
# Handle the degenerate case
40+
if iszero(k)
41+
return float(last(promote(x, x < 0)))
42+
end
43+
return last(gamma_inc(k, max(0, x)/θ))
44+
end
45+
gammaccdf(k::Real, θ::Real, x::Real) = gammaccdf(map(float, promote(k, θ, x))...)
46+
47+
gammalogcdf(k::Real, θ::Real, x::Real) = _gammalogcdf(map(float, promote(k, θ, x))...)
48+
49+
# Implemented via the non-log version. For tiny values, we recompute the result with
50+
# loggamma. In that situation, there is little risk of significant cancellation.
51+
function _gammalogcdf(k::Float64, θ::Float64, x::Float64)
52+
# Handle the degenerate case
53+
if iszero(k)
54+
return log(x >= 0)
55+
end
56+
57+
xdθ = max(0, x)/θ
58+
l, u = gamma_inc(k, xdθ)
59+
if l < floatmin(Float64) && isfinite(k) && isfinite(xdθ)
60+
return -log(k) + k*log(xdθ) - xdθ + log(drummond1F1(1.0, 1 + k, xdθ)) - loggamma(k)
61+
elseif l < 0.7
62+
return log(l)
63+
else
64+
return log1p(-u)
65+
end
66+
end
67+
function _gammalogcdf(k::T, θ::T, x::T) where {T<:Union{Float16,Float32}}
68+
return T(_gammalogcdf(Float64(k), Float64(θ), Float64(x)))
69+
end
70+
71+
gammalogccdf(k::Real, θ::Real, x::Real) = _gammalogccdf(map(float, promote(k, θ, x))...)
72+
73+
# Implemented via the non-log version. For tiny values, we recompute the result with
74+
# loggamma. In that situation, there is little risk of significant cancellation.
75+
function _gammalogccdf(k::Float64, θ::Float64, x::Float64)
76+
# Handle the degenerate case
77+
if iszero(k)
78+
return log(x < 0)
79+
end
80+
81+
xdθ = max(0, x)/θ
82+
l, u = gamma_inc(k, xdθ)
83+
if u < floatmin(Float64)
84+
return loggamma(k, xdθ) - loggamma(k)
85+
elseif u < 0.7
86+
return log(u)
87+
else
88+
return log1p(-l)
89+
end
90+
end
91+
92+
function _gammalogccdf(k::T, θ::T, x::T) where {T<:Union{Float16,Float32}}
93+
return T(_gammalogccdf(Float64(k), Float64(θ), Float64(x)))
94+
end
95+
96+
function gammainvcdf(k::Real, θ::Real, p::Real)
97+
_k, _θ, _p = promote(k, θ, p)
98+
return _θ*gamma_inc_inv(_k, _p, 1 - _p)
99+
end
100+
101+
function gammainvccdf(k::Real, θ::Real, p::Real)
102+
_k, _θ, _p = promote(k, θ, p)
103+
return _θ*gamma_inc_inv(_k, 1 - _p, _p)
104+
end

‎src/distrs/pois.jl‎

Lines changed: 15 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,13 @@
11
# functions related to Poisson distribution
22

33
# R implementations
4-
# For pdf and logpdf we use the Julia implementation
54
using .RFunctions:
6-
poiscdf,
7-
poisccdf,
8-
poislogcdf,
9-
poislogccdf,
5+
# poispdf,
6+
# poislogpdf,
7+
# poiscdf,
8+
# poisccdf,
9+
# poislogcdf,
10+
# poislogccdf,
1011
poisinvcdf,
1112
poisinvccdf,
1213
poisinvlogcdf,
@@ -20,3 +21,12 @@ function poislogpdf(λ::T, x::T) where {T <: Real}
2021
val = xlogy(x, λ) - λ - loggamma(x + 1)
2122
return x >= 0 && isinteger(x) ? val : oftype(val, -Inf)
2223
end
24+
25+
# Just use the Gamma definitions
26+
poiscdf(λ::Real, x::Real) = gammaccdf(max(0, floor(x + 1)), 1, λ)
27+
28+
poisccdf(λ::Real, x::Real) = gammacdf(max(0, floor(x + 1)), 1, λ)
29+
30+
poislogcdf(λ::Real, x::Real) = gammalogccdf(max(0, floor(x + 1)), 1, λ)
31+
32+
poislogccdf(λ::Real, x::Real) = gammalogcdf(max(0, floor(x + 1)), 1, λ)

‎src/distrs/tdist.jl‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,9 @@
11
# functions related to student's T distribution
22

33
# R implementations
4-
# For pdf and logpdf we use the Julia implementation
54
using .RFunctions:
5+
# tdistpdf,
6+
# tdistlogpdf,
67
tdistcdf,
78
tdistccdf,
89
tdistlogcdf,

0 commit comments

Comments
 (0)