WIP: remove Rmath dependence

src/Distributions.jl

@@ -217,7 +217,10 @@ include("constants.jl")
 
 include("fallbacks.jl")
 include("rmath.jl")
-include("specialfuns.jl")
+include("specialfuns/misc.jl")
+include("specialfuns/normal.jl")
+include("specialfuns/log.jl")
+include("specialfuns/gammabeta.jl")
 include("tvpack.jl")
 include("utils.jl")
 
@@ -229,6 +232,7 @@ include(joinpath("univariate", "bernoulli.jl"))
 include(joinpath("univariate", "beta.jl"))
 include(joinpath("univariate", "betaprime.jl"))
 include(joinpath("univariate", "binomial.jl"))
+include(joinpath("univariate", "binomial_rand.jl"))
 include(joinpath("univariate", "categorical.jl"))
 include(joinpath("univariate", "cauchy.jl"))
 include(joinpath("univariate", "chi.jl"))
@@ -262,6 +266,7 @@ include(joinpath("univariate", "normal.jl"))
 include(joinpath("univariate", "normalcanon.jl"))
 include(joinpath("univariate", "pareto.jl"))
 include(joinpath("univariate", "poisson.jl"))
+include(joinpath("univariate", "poisson_rand.jl"))
 include(joinpath("univariate", "rayleigh.jl"))
 include(joinpath("univariate", "skellam.jl"))
 include(joinpath("univariate", "tdist.jl"))

src/constants.jl

@@ -3,6 +3,9 @@
 import Base.@math_const
 
 @math_const twoπ   6.2831853071795864769 big(2.) * π
+@math_const halfπ  1.5707963267948966192 big(0.5) * π
 @math_const √2     1.4142135623730950488 sqrt(big(2.))
 @math_const log2π  1.8378770664093454836 log(big(2.)*π)
+@math_const hlog2π 0.9189385332046727418 0.5*log(big(2.)*π)
 @math_const √2π    2.5066282746310005024 sqrt(big(2.)*π)
+@math_const r√2π   0.3989422804014326779 1/sqrt(big(2.)*π)

src/fallbacks.jl

@@ -105,6 +105,47 @@ logccdf(d::Distribution, q::Real) = log(ccdf(d,q))
 invlogccdf(d::Distribution, lp::Real) = quantile(d, -expm1(lp))
 invlogcdf(d::Distribution, lp::Real) = quantile(d, exp(lp))
 
+function quantile(d::ContinuousUnivariateDistribution, α::Real)
+
+    if α < 0 || α > 1 return NaN end
+    if α == 0 return 0.0 end
+    if α == 1 return 1.0 end
+
+    cc = 10eps()
+    e = sqrt(eps())
+    x = mode(d)
+    while true
+        dx = (cdf(d, x)::Float64 - α)/max(e,pdf(d, x)::Float64)
+        if abs(dx) < cc 
+            x -= dx
+            return x
+        end
+        t = x - dx
+        while !insupport(d, t)
+            dx *= 0.5
+            t = x - dx
+        end
+        x = t
+    end
+end
+
+function quantile(d::DiscreteUnivariateDistribution, α::Real)
+    if α < 0 || α > 1 return NaN end
+    if α == 0 return 0 end
+    if α == 1 return d.size end
+    qc = itrunc(quantile(Normal(mean(d), std(d)), α))
+    if α < cdf(d, qc)
+        qc -= 1
+        while α < cdf(d, qc)
+            qc -= 1
+        end
+        return qc + 1
+    end
+    while α > cdf(d, qc)
+        qc += 1
+    end
+    return qc
+end
 
 #### log likelihood ####