Remove Rmath sampling methods

Project.toml

@@ -1,6 +1,6 @@
 name = "StatsFuns"
 uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
-version = "1.5.2"
+version = "2.0.0"
 
 [deps]
 HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"

src/StatsFuns.jl

@@ -2,6 +2,7 @@ module StatsFuns
 
 using Base: Math.@horner
 using Reexport: @reexport
+using Rmath: Rmath
 using SpecialFunctions: beta_inc, beta_inc_inv, digamma,
     erfc, erfcinv, erfcx, gamma_inc, gamma_inc_inv, logbeta, loggamma
 
@@ -264,11 +265,8 @@ export
 
 ## source files
 include("misc.jl")
-include("rmath.jl")
 include("tvpack.jl")
 
-using .RFunctions
-
 include("distrs/beta.jl")
 include("distrs/binom.jl")
 include("distrs/chisq.jl")

src/distrs/beta.jl

@@ -2,19 +2,6 @@
 
 using HypergeometricFunctions: _₂F₁
 
-# R implementations
-using .RFunctions:
-    # betapdf,
-    # betalogpdf,
-    # betacdf,
-    # betaccdf,
-    # betalogcdf,
-    # betalogccdf,
-    # betainvcdf,
-    # betainvccdf,
-    betainvlogcdf,
-    betainvlogccdf
-
 # Julia implementations
 betapdf(α::Real, β::Real, x::Real) = exp(betalogpdf(α, β, x))
 
@@ -113,3 +100,13 @@ function betainvccdf(α::Real, β::Real, p::Real)
 
     return last(beta_inc_inv(β, α, p))
 end
+
+# Rmath implementations
+function betainvlogcdf(α::Real, β::Real, lq::Real)
+    T = float(Base.promote_typeof(α, β, lq))
+    return convert(T, Rmath.qbeta(lq, α, β, true, true))
+end
+function betainvlogccdf(α::Real, β::Real, lq::Real)
+    T = float(Base.promote_typeof(α, β, lq))
+    return convert(T, Rmath.qbeta(lq, α, β, false, true))
+end

src/distrs/binom.jl

@@ -1,18 +1,5 @@
 # functions related to binomial distribution
 
-# R implementations
-using .RFunctions:
-    # binompdf,
-    # binomlogpdf,
-    # binomcdf,
-    # binomccdf,
-    # binomlogcdf,
-    # binomlogccdf,
-    binominvcdf,
-    binominvccdf,
-    binominvlogcdf,
-    binominvlogccdf
-
 # Julia implementations
 binompdf(n::Real, p::Real, k::Real) = exp(binomlogpdf(n, p, k))
 
@@ -42,3 +29,21 @@ for l in ("", "log"), compl in (false, true)
         end
     end
 end
+
+# Rmath implementations
+function binominvcdf(n::Real, p::Real, q::Real)
+    T = float(Base.promote_typeof(n, p, q))
+    return convert(T, Rmath.qbinom(q, n, p, true, false))
+end
+function binominvccdf(n::Real, p::Real, q::Real)
+    T = float(Base.promote_typeof(n, p, q))
+    return convert(T, Rmath.qbinom(q, n, p, false, false))
+end
+function binominvlogcdf(n::Real, p::Real, lq::Real)
+    T = float(Base.promote_typeof(n, p, lq))
+    return convert(T, Rmath.qbinom(lq, n, p, true, true))
+end
+function binominvlogccdf(n::Real, p::Real, lq::Real)
+    T = float(Base.promote_typeof(n, p, lq))
+    return convert(T, Rmath.qbinom(lq, n, p, false, true))
+end

src/distrs/gamma.jl

@@ -2,19 +2,6 @@
 
 using HypergeometricFunctions: _₁F₁
 
-# R implementations
-using .RFunctions:
-    # gammapdf,
-    # gammalogpdf,
-    # gammacdf,
-    # gammaccdf,
-    # gammalogcdf,
-    # gammalogccdf,
-    # gammainvcdf,
-    # gammainvccdf,
-    gammainvlogcdf,
-    gammainvlogccdf
-
 # Julia implementations
 gammapdf(k::Real, θ::Real, x::Real) = exp(gammalogpdf(k, θ, x))
 
@@ -107,3 +94,13 @@ function gammainvccdf(k::Real, θ::Real, p::Real)
     _k, _θ, _p = promote(k, θ, p)
     return _θ * gamma_inc_inv(_k, 1 - _p, _p)
 end
+
+# Rmath implementations
+function gammainvlogcdf(k::Real, θ::Real, lq::Real)
+    T = float(Base.promote_typeof(k, θ, lq))
+    return convert(T, Rmath.qgamma(lq, k, θ, true, true))
+end
+function gammainvlogccdf(k::Real, θ::Real, lq::Real)
+    T = float(Base.promote_typeof(k, θ, lq))
+    return convert(T, Rmath.qgamma(lq, k, θ, false, true))
+end

src/distrs/hyper.jl

@@ -1,14 +1,45 @@
 # functions related to hyper-geometric distribution
 
-# R implementations
-using .RFunctions:
-    hyperpdf,
-    hyperlogpdf,
-    hypercdf,
-    hyperccdf,
-    hyperlogcdf,
-    hyperlogccdf,
-    hyperinvcdf,
-    hyperinvccdf,
-    hyperinvlogcdf,
-    hyperinvlogccdf
+# Rmath implementations
+function hyperpdf(ms::Real, mf::Real, n::Real, x::Real)
+    T = float(Base.promote_typeof(ms, mf, n, x))
+    return convert(T, Rmath.dhyper(x, ms, mf, n, false))
+end
+function hyperlogpdf(ms::Real, mf::Real, n::Real, x::Real)
+    T = float(Base.promote_typeof(ms, mf, n, x))
+    return convert(T, Rmath.dhyper(x, ms, mf, n, true))
+end
+
+function hypercdf(ms::Real, mf::Real, n::Real, x::Real)
+    T = float(Base.promote_typeof(ms, mf, n, x))
+    return convert(T, Rmath.phyper(x, ms, mf, n, true, false))
+end
+function hyperccdf(ms::Real, mf::Real, n::Real, x::Real)
+    T = float(Base.promote_typeof(ms, mf, n, x))
+    return convert(T, Rmath.phyper(x, ms, mf, n, false, false))
+end
+function hyperlogcdf(ms::Real, mf::Real, n::Real, x::Real)
+    T = float(Base.promote_typeof(ms, mf, n, x))
+    return convert(T, Rmath.phyper(x, ms, mf, n, true, true))
+end
+function hyperlogccdf(ms::Real, mf::Real, n::Real, x::Real)
+    T = float(Base.promote_typeof(ms, mf, n, x))
+    return convert(T, Rmath.phyper(x, ms, mf, n, false, true))
+end
+
+function hyperinvcdf(ms::Real, mf::Real, n::Real, q::Real)
+    T = float(Base.promote_typeof(ms, mf, n, q))
+    return convert(T, Rmath.qhyper(q, ms, mf, n, true, false))
+end
+function hyperinvccdf(ms::Real, mf::Real, n::Real, q::Real)
+    T = float(Base.promote_typeof(ms, mf, n, q))
+    return convert(T, Rmath.qhyper(q, ms, mf, n, false, false))
+end
+function hyperinvlogcdf(ms::Real, mf::Real, n::Real, lq::Real)
+    T = float(Base.promote_typeof(ms, mf, n, lq))
+    return convert(T, Rmath.qhyper(lq, ms, mf, n, true, true))
+end
+function hyperinvlogccdf(ms::Real, mf::Real, n::Real, lq::Real)
+    T = float(Base.promote_typeof(ms, mf, n, lq))
+    return convert(T, Rmath.qhyper(lq, ms, mf, n, false, true))
+end

src/distrs/nbeta.jl

@@ -1,14 +1,45 @@
 # functions related to noncentral beta distribution
 
-# R implementations
-using .RFunctions:
-    nbetapdf,
-    nbetalogpdf,
-    nbetacdf,
-    nbetaccdf,
-    nbetalogcdf,
-    nbetalogccdf,
-    nbetainvcdf,
-    nbetainvccdf,
-    nbetainvlogcdf,
-    nbetainvlogccdf
+# Rmath implementations
+function nbetapdf(α::Real, β::Real, λ::Real, x::Real)
+    T = float(Base.promote_typeof(α, β, λ, x))
+    return convert(T, Rmath.dnbeta(x, α, β, λ, false))
+end
+function nbetalogpdf(α::Real, β::Real, λ::Real, x::Real)
+    T = float(Base.promote_typeof(α, β, λ, x))
+    return convert(T, Rmath.dnbeta(x, α, β, λ, true))
+end
+
+function nbetacdf(α::Real, β::Real, λ::Real, x::Real)
+    T = float(Base.promote_typeof(α, β, λ, x))
+    return convert(T, Rmath.pnbeta(x, α, β, λ, true, false))
+end
+function nbetaccdf(α::Real, β::Real, λ::Real, x::Real)
+    T = float(Base.promote_typeof(α, β, λ, x))
+    return convert(T, Rmath.pnbeta(x, α, β, λ, false, false))
+end
+function nbetalogcdf(α::Real, β::Real, λ::Real, x::Real)
+    T = float(Base.promote_typeof(α, β, λ, x))
+    return convert(T, Rmath.pnbeta(x, α, β, λ, true, true))
+end
+function nbetalogccdf(α::Real, β::Real, λ::Real, x::Real)
+    T = float(Base.promote_typeof(α, β, λ, x))
+    return convert(T, Rmath.pnbeta(x, α, β, λ, false, true))
+end
+
+function nbetainvcdf(α::Real, β::Real, λ::Real, q::Real)
+    T = float(Base.promote_typeof(α, β, λ, q))
+    return convert(T, Rmath.qnbeta(q, α, β, λ, true, false))
+end
+function nbetainvccdf(α::Real, β::Real, λ::Real, q::Real)
+    T = float(Base.promote_typeof(α, β, λ, q))
+    return convert(T, Rmath.qnbeta(q, α, β, λ, false, false))
+end
+function nbetainvlogcdf(α::Real, β::Real, λ::Real, lq::Real)
+    T = float(Base.promote_typeof(α, β, λ, lq))
+    return convert(T, Rmath.qnbeta(lq, α, β, λ, true, true))
+end
+function nbetainvlogccdf(α::Real, β::Real, λ::Real, lq::Real)
+    T = float(Base.promote_typeof(α, β, λ, lq))
+    return convert(T, Rmath.qnbeta(lq, α, β, λ, false, true))
+end

src/distrs/nbinom.jl

@@ -1,14 +1,45 @@
 # functions related to negative binomial distribution
 
-# R implementations
-using .RFunctions:
-    nbinompdf,
-    nbinomlogpdf,
-    nbinomcdf,
-    nbinomccdf,
-    nbinomlogcdf,
-    nbinomlogccdf,
-    nbinominvcdf,
-    nbinominvccdf,
-    nbinominvlogcdf,
-    nbinominvlogccdf
+# Rmath implementations
+function nbinompdf(r::Real, p::Real, x::Real)
+    T = float(Base.promote_typeof(r, p, x))
+    return convert(T, Rmath.dnbinom(x, r, p, false))
+end
+function nbinomlogpdf(r::Real, p::Real, x::Real)
+    T = float(Base.promote_typeof(r, p, x))
+    return convert(T, Rmath.dnbinom(x, r, p, true))
+end
+
+function nbinomcdf(r::Real, p::Real, x::Real)
+    T = float(Base.promote_typeof(r, p, x))
+    return convert(T, Rmath.pnbinom(x, r, p, true, false))
+end
+function nbinomccdf(r::Real, p::Real, x::Real)
+    T = float(Base.promote_typeof(r, p, x))
+    return convert(T, Rmath.pnbinom(x, r, p, false, false))
+end
+function nbinomlogcdf(r::Real, p::Real, x::Real)
+    T = float(Base.promote_typeof(r, p, x))
+    return convert(T, Rmath.pnbinom(x, r, p, true, true))
+end
+function nbinomlogccdf(r::Real, p::Real, x::Real)
+    T = float(Base.promote_typeof(r, p, x))
+    return convert(T, Rmath.pnbinom(x, r, p, false, true))
+end
+
+function nbinominvcdf(r::Real, p::Real, q::Real)
+    T = float(Base.promote_typeof(r, p, q))
+    return convert(T, Rmath.qnbinom(q, r, p, true, false))
+end
+function nbinominvccdf(r::Real, p::Real, q::Real)
+    T = float(Base.promote_typeof(r, p, q))
+    return convert(T, Rmath.qnbinom(q, r, p, false, false))
+end
+function nbinominvlogcdf(r::Real, p::Real, lq::Real)
+    T = float(Base.promote_typeof(r, p, lq))
+    return convert(T, Rmath.qnbinom(lq, r, p, true, true))
+end
+function nbinominvlogccdf(r::Real, p::Real, lq::Real)
+    T = float(Base.promote_typeof(r, p, lq))
+    return convert(T, Rmath.qnbinom(lq, r, p, false, true))
+end

src/distrs/nchisq.jl

@@ -1,14 +1,45 @@
 # functions related to noncentral chi-square distribution
 
-# R implementations
-using .RFunctions:
-    nchisqpdf,
-    nchisqlogpdf,
-    nchisqcdf,
-    nchisqccdf,
-    nchisqlogcdf,
-    nchisqlogccdf,
-    nchisqinvcdf,
-    nchisqinvccdf,
-    nchisqinvlogcdf,
-    nchisqinvlogccdf
+# Rmath implementations
+function nchisqpdf(k::Real, λ::Real, x::Real)
+    T = float(Base.promote_typeof(k, λ, x))
+    return convert(T, Rmath.dnchisq(x, k, λ, false))
+end
+function nchisqlogpdf(k::Real, λ::Real, x::Real)
+    T = float(Base.promote_typeof(k, λ, x))
+    return convert(T, Rmath.dnchisq(x, k, λ, true))
+end
+
+function nchisqcdf(k::Real, λ::Real, x::Real)
+    T = float(Base.promote_typeof(k, λ, x))
+    return convert(T, Rmath.pnchisq(x, k, λ, true, false))
+end
+function nchisqccdf(k::Real, λ::Real, x::Real)
+    T = float(Base.promote_typeof(k, λ, x))
+    return convert(T, Rmath.pnchisq(x, k, λ, false, false))
+end
+function nchisqlogcdf(k::Real, λ::Real, x::Real)
+    T = float(Base.promote_typeof(k, λ, x))
+    return convert(T, Rmath.pnchisq(x, k, λ, true, true))
+end
+function nchisqlogccdf(k::Real, λ::Real, x::Real)
+    T = float(Base.promote_typeof(k, λ, x))
+    return convert(T, Rmath.pnchisq(x, k, λ, false, true))
+end
+
+function nchisqinvcdf(k::Real, λ::Real, q::Real)
+    T = float(Base.promote_typeof(k, λ, q))
+    return convert(T, Rmath.qnchisq(q, k, λ, true, false))
+end
+function nchisqinvccdf(k::Real, λ::Real, q::Real)
+    T = float(Base.promote_typeof(k, λ, q))
+    return convert(T, Rmath.qnchisq(q, k, λ, false, false))
+end
+function nchisqinvlogcdf(k::Real, λ::Real, lq::Real)
+    T = float(Base.promote_typeof(k, λ, lq))
+    return convert(T, Rmath.qnchisq(lq, k, λ, true, true))
+end
+function nchisqinvlogccdf(k::Real, λ::Real, lq::Real)
+    T = float(Base.promote_typeof(k, λ, lq))
+    return convert(T, Rmath.qnchisq(lq, k, λ, false, true))
+end