diff --git a/.appveyor.yml b/.appveyor.yml index deb2424..da92a6a 100644 --- a/.appveyor.yml +++ b/.appveyor.yml @@ -1,7 +1,7 @@ environment: matrix: - - julia_version: 0.7 - julia_version: 1 + - julia_version: 1.2 - julia_version: nightly platform: diff --git a/.travis.yml b/.travis.yml index f366abf..4e88b8e 100644 --- a/.travis.yml +++ b/.travis.yml @@ -3,8 +3,8 @@ os: - linux - osx julia: - - 0.7 - 1.0 + - 1.2 - nightly notifications: email: false diff --git a/Project.toml b/Project.toml new file mode 100644 index 0000000..30888fc --- /dev/null +++ b/Project.toml @@ -0,0 +1,19 @@ +name = "StatsFuns" +uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c" +version = "0.9.0" + +[deps] +Rmath = "79098fc4-a85e-5d69-aa6a-4863f24498fa" +SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" + +[compat] +Rmath = "0.4, 0.5" +SpecialFunctions = "0.8" +julia = "1" + +[extras] +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[targets] +test = ["ForwardDiff", "Test"] diff --git a/REQUIRE b/REQUIRE deleted file mode 100644 index 02ca4f6..0000000 --- a/REQUIRE +++ /dev/null @@ -1,3 +0,0 @@ -julia 0.7-beta -Rmath 0.4.0 -SpecialFunctions diff --git a/src/distrs/beta.jl b/src/distrs/beta.jl index 2a8919a..3d1bfae 100644 --- a/src/distrs/beta.jl +++ b/src/distrs/beta.jl @@ -16,4 +16,4 @@ import .RFunctions: betapdf(α::Real, β::Real, x::Number) = x^(α - 1) * (1 - x)^(β - 1) / beta(α, β) # logpdf for numbers with generic types -betalogpdf(α::Real, β::Real, x::Number) = (α - 1) * log(x) + (β - 1) * log1p(-x) - lbeta(α, β) +betalogpdf(α::Real, β::Real, x::Number) = (α - 1) * log(x) + (β - 1) * log1p(-x) - logbeta(α, β) diff --git a/src/distrs/binom.jl b/src/distrs/binom.jl index 92b440c..e1237f1 100644 --- a/src/distrs/binom.jl +++ b/src/distrs/binom.jl @@ -16,4 +16,4 @@ import .RFunctions: binompdf(n::Real, p::Real, k::Real) = exp(binomlogpdf(n, p, k)) # logpdf for numbers with generic types -binomlogpdf(n::Real, p::Real, k::Real) = -log1p(n) - lbeta(n - k + 1, k + 1) + k * log(p) + (n - k) * log1p(-p) +binomlogpdf(n::Real, p::Real, k::Real) = -log1p(n) - logbeta(n - k + 1, k + 1) + k * log(p) + (n - k) * log1p(-p) diff --git a/src/distrs/chisq.jl b/src/distrs/chisq.jl index 17e4987..606923e 100644 --- a/src/distrs/chisq.jl +++ b/src/distrs/chisq.jl @@ -21,5 +21,5 @@ end # logpdf for numbers with generic types function chisqlogpdf(k::Real, x::Number) hk = k / 2 # half k - -hk * log(oftype(hk, 2)) - lgamma(hk) + (hk - 1) * log(x) - x / 2 + -hk * log(oftype(hk, 2)) - loggamma(hk) + (hk - 1) * log(x) - x / 2 end diff --git a/src/distrs/fdist.jl b/src/distrs/fdist.jl index b46c796..b3041e3 100644 --- a/src/distrs/fdist.jl +++ b/src/distrs/fdist.jl @@ -16,4 +16,4 @@ import .RFunctions: fdistpdf(ν1::Real, ν2::Real, x::Number) = sqrt((ν1 * x)^ν1 * ν2^ν2 / (ν1 * x + ν2)^(ν1 + ν2)) / (x * beta(ν1 / 2, ν2 / 2)) # logpdf for numbers with generic types -fdistlogpdf(ν1::Real, ν2::Real, x::Number) = (ν1 * log(ν1 * x) + ν2 * log(ν2) - (ν1 + ν2) * log(ν1 * x + ν2)) / 2 - log(x) - lbeta(ν1 / 2, ν2 / 2) +fdistlogpdf(ν1::Real, ν2::Real, x::Number) = (ν1 * log(ν1 * x) + ν2 * log(ν2) - (ν1 + ν2) * log(ν1 * x + ν2)) / 2 - log(x) - logbeta(ν1 / 2, ν2 / 2) diff --git a/src/distrs/gamma.jl b/src/distrs/gamma.jl index 0ad808c..26dbef6 100644 --- a/src/distrs/gamma.jl +++ b/src/distrs/gamma.jl @@ -16,4 +16,4 @@ import .RFunctions: gammapdf(k::Real, θ::Real, x::Number) = 1 / (gamma(k) * θ^k) * x^(k - 1) * exp(-x / θ) # logpdf for numbers with generic types -gammalogpdf(k::Real, θ::Real, x::Number) = -lgamma(k) - k * log(θ) + (k - 1) * log(x) - x / θ +gammalogpdf(k::Real, θ::Real, x::Number) = -loggamma(k) - k * log(θ) + (k - 1) * log(x) - x / θ diff --git a/src/distrs/pois.jl b/src/distrs/pois.jl index adedaa6..1c1ad3b 100644 --- a/src/distrs/pois.jl +++ b/src/distrs/pois.jl @@ -15,7 +15,7 @@ import .RFunctions: # generic versions poispdf(λ::Real, x::Real) = exp(poislogpdf(λ, x)) -poislogpdf(λ::T, x::T) where {T <: Real} = xlogy(x, λ) - λ - lgamma(x + 1) +poislogpdf(λ::T, x::T) where {T <: Real} = xlogy(x, λ) - λ - loggamma(x + 1) poislogpdf(λ::Number, x::Number) = poislogpdf(promote(float(λ), x)...) diff --git a/src/distrs/tdist.jl b/src/distrs/tdist.jl index 21f1fa3..322a103 100644 --- a/src/distrs/tdist.jl +++ b/src/distrs/tdist.jl @@ -16,4 +16,4 @@ import .RFunctions: tdistpdf(ν::Real, x::Number) = gamma((ν + 1) / 2) / (sqrt(ν * pi) * gamma(ν / 2)) * (1 + x^2 / ν)^(-(ν + 1) / 2) # logpdf for numbers with generic types -tdistlogpdf(ν::Real, x::Number) = lgamma((ν + 1) / 2) - log(ν * pi) / 2 - lgamma(ν / 2) + (-(ν + 1) / 2) * log(1 + x^2 / ν) +tdistlogpdf(ν::Real, x::Number) = loggamma((ν + 1) / 2) - log(ν * pi) / 2 - loggamma(ν / 2) + (-(ν + 1) / 2) * log(1 + x^2 / ν) diff --git a/src/misc.jl b/src/misc.jl index 2177843..8e76c2e 100644 --- a/src/misc.jl +++ b/src/misc.jl @@ -9,7 +9,7 @@ function logmvgamma(p::Int, a::Real) # NOTE: one(a) factors are here to prevent unnecessary promotion of Float32 res = p * (p - 1) * log(pi * one(a)) / 4 for ii in 1:p - res += lgamma(a + (1 - ii) * one(a)/ 2) + res += loggamma(a + (1 - ii) * one(a) / 2) end return res end @@ -19,7 +19,7 @@ end The remainder term after [Stirling's approximation](https://en.wikipedia.org/wiki/Stirling%27s_approximation) -to [`lgamma`](@ref): +to [`loggamma`](@ref): ```math \\log \\Gamma(x) \\approx x \\log(x) - x + \\log(2\\pi/x)/2 = \\log(x)*(x-1/2) + \\log(2\\pi)/2 - x @@ -27,7 +27,7 @@ to [`lgamma`](@ref): In Julia syntax, this means: - lstirling_asym(x) = lgamma(x) + x - (x-0.5)*log(x) - 0.5*log(2π) + lstirling_asym(x) = loggamma(x) + x - (x-0.5)*log(x) - 0.5*log(2π) For sufficiently large `x`, this can be approximated using the asymptotic _Stirling's series_ ([DLMF 5.11.1](https://dlmf.nist.gov/5.11.1)): @@ -53,7 +53,7 @@ which is < 1/2 ulp for x >= 10.0, and total numeric error appears to be < 2 ulps """ function lstirling_asym end -lstirling_asym(x::BigFloat) = lgamma(x) + x - log(x)*(x - big(0.5)) - log2π/big(2) +lstirling_asym(x::BigFloat) = loggamma(x) + x - log(x)*(x - big(0.5)) - log2π/big(2) lstirling_asym(x::Integer) = lstirling_asym(float(x)) diff --git a/test/REQUIRE b/test/REQUIRE deleted file mode 100644 index ea9c38c..0000000 --- a/test/REQUIRE +++ /dev/null @@ -1 +0,0 @@ -ForwardDiff 0.3.3 diff --git a/test/misc.jl b/test/misc.jl index ab4c49e..c850d17 100644 --- a/test/misc.jl +++ b/test/misc.jl @@ -4,22 +4,23 @@ using SpecialFunctions, StatsFuns @testset "type behavior" for eltya in (Float32, Float64) p = rand(1:50) a = rand(eltya) - @test typeof(logmvgamma(p, a)) == eltya + # add p since loggamma is only define for positive arguments + @test typeof(logmvgamma(p, a + p)) == eltya end - @testset "consistent with lgamma" for eltya in (Float32, Float64) + @testset "consistent with loggamma" for eltya in (Float32, Float64) # Γ₁(a) = Γ(a), Γ₂(a) = √π Γ(a) Γ(a - 0.5), etc - a = rand(eltya) - @test logmvgamma(1, a) ≈ lgamma(a) - @test logmvgamma(2, a) ≈ eltya(0.5logπ) + lgamma(a) + lgamma(a - eltya(0.5)) - @test logmvgamma(3, a) ≈ eltya((3/2)*logπ) + lgamma(a) + lgamma(a - eltya(0.5)) + lgamma(a - one(a)) + a = rand(eltya) + 1 # add one since loggamma is only define for positive arguments + @test logmvgamma(1, a) ≈ loggamma(a) + @test logmvgamma(2, a) ≈ eltya(0.5logπ) + loggamma(a) + loggamma(a - eltya(0.5)) + @test logmvgamma(3, a) ≈ eltya((3/2)*logπ) + loggamma(a) + loggamma(a - eltya(0.5)) + loggamma(a - one(a)) end @testset "consistent with itself" for eltya in (Float32, Float64) # Γᵢ(a) = (π^{i-1/2}) Γ(a) Γᵢ₋₁(a - 0.5) for p in 1:50 - a = rand(eltya) - @test logmvgamma(p, a) ≈ eltya((p/2-1/2)*logπ) + lgamma(a) + logmvgamma(p - 1, a - eltya(0.5)) + a = rand(eltya) + p # add p since loggamma is only define for positive arguments + @test logmvgamma(p, a) ≈ eltya((p/2-1/2)*logπ) + loggamma(a) + logmvgamma(p - 1, a - eltya(0.5)) end end end