Merge branch 'master' into dw/pairwise

Author: devmotion

Project.toml

@@ -1,7 +1,7 @@
 name = "Distributions"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
 authors = ["JuliaStats"]
-version = "0.25.28"
+version = "0.25.29"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

docs/src/univariate.md

@@ -327,6 +327,13 @@ LogNormal
 plotdensity((0, 5), LogNormal, (0, 1)) # hide
 ```
 
+```@docs
+LogUniform
+```
+```@example plotdensity
+plotdensity((0, 11), LogUniform, (1, 10)) # hide
+```
+
 ```@docs
 NoncentralBeta
 ```

src/Distributions.jl

@@ -118,6 +118,7 @@ export
     LocationScale,
     Logistic,
     LogNormal,
+    LogUniform,
     LogitNormal,
     MatrixBeta,
     MatrixFDist,

src/truncate.jl

@@ -172,3 +172,4 @@ _use_multline_show(d::Truncated) = _use_multline_show(d.untruncated)
 include(joinpath("truncated", "normal.jl"))
 include(joinpath("truncated", "exponential.jl"))
 include(joinpath("truncated", "uniform.jl"))
+include(joinpath("truncated", "loguniform.jl"))

src/truncated/loguniform.jl

@@ -0,0 +1 @@
+truncated(d::LogUniform, lo::T, hi::T) where {T<:Real} = LogUniform(max(d.a, lo), min(d.b, hi))

src/univariate/continuous/loguniform.jl

@@ -0,0 +1,75 @@
+"""
+    LogUniform(a,b)
+
+A positive random variable `X` is log-uniformly with parameters `a` and `b` if the logarithm of `X` is `Uniform(log(a), log(b))`.
+The *log uniform* distribution is also known as *reciprocal distribution*.
+```julia
+LogUniform(1,10)
+```
+External links
+
+* [Log uniform distribution on Wikipedia](https://en.wikipedia.org/wiki/Reciprocal_distribution)
+"""
+struct LogUniform{T<:Real} <: ContinuousUnivariateDistribution
+    a::T
+    b::T
+    LogUniform{T}(a::T, b::T) where {T <: Real} = new{T}(a, b)
+end
+
+function LogUniform(a::T, b::T; check_args=true) where {T <: Real}
+    check_args && @check_args(LogUniform, 0 < a < b)
+    LogUniform{T}(a, b)
+end
+
+LogUniform(a::Real, b::Real; kwargs...) = LogUniform(promote(a, b)...; kwargs...)
+
+convert(::Type{LogUniform{T}}, d::LogUniform) where {T<:Real} = LogUniform(T(d.a), T(d.b))
+Base.minimum(d::LogUniform) = d.a
+Base.maximum(d::LogUniform) = d.b
+
+#### Parameters
+params(d::LogUniform) = (d.a, d.b)
+partype(::LogUniform{T}) where {T<:Real} = T
+
+#### Statistics
+
+function mean(d::LogUniform)
+    a, b = params(d)
+    (b - a) / log(b/a)
+end
+function var(d::LogUniform)
+    a, b = params(d)
+    log_ba = log(b/a)
+    (b^2 - a^2) / (2*log_ba) - ((b-a)/ log_ba)^2
+end
+mode(d::LogUniform)   = d.a
+modes(d::LogUniform)  = partype(d)[]
+
+function entropy(d::LogUniform)
+    a,b = params(d)
+    log(a * b) / 2 + log(log(b / a))
+end
+#### Evaluation
+function pdf(d::LogUniform, x::Real)
+    x1, a, b = promote(x, params(d)...) # ensure e.g. pdf(LogUniform(1,2), 1f0)::Float32
+    res = inv(x1 * log(b / a))
+    return insupport(d, x1) ? res : zero(res)
+end
+function cdf(d::LogUniform, x::Real)
+    x1, a, b = promote(x, params(d)...) # ensure e.g. cdf(LogUniform(1,2), 1f0)::Float32
+    x1 = clamp(x1, a, b)
+    return log(x1 / a) / log(b / a)
+end
+logpdf(d::LogUniform, x::Real) = log(pdf(d,x))
+
+function quantile(d::LogUniform, p::Real)
+    p1,a,b = promote(p, params(d)...) # ensure e.g. quantile(LogUniform(1,2), 1f0)::Float32
+    exp(p1 * log(b/a)) * a
+end
+
+function kldivergence(p::LogUniform, q::LogUniform)
+    ap, bp, aq, bq = promote(params(p)..., params(q)...)
+    finite = aq <= ap < bp <= bq
+    res = log(log(bq / aq) / log(bp / ap))
+    return finite ? res : oftype(res, Inf)
+end

src/univariates.jl

@@ -716,6 +716,7 @@ const continuous_distributions = [
     "triangular",
     "triweight",
     "uniform",
+    "loguniform", # depends on Uniform
     "vonmises",
     "weibull"
 ]

test/loguniform.jl

@@ -0,0 +1,85 @@
+module TestLogUniform
+using Test
+using Distributions
+import Random
+
+@testset "LogUniform" begin
+    rng = Random.MersenneTwister(0)
+
+    @test pdf(LogUniform(1f0, 2f0), 1) isa Float32
+    @test pdf(LogUniform(1, 2), 1f0) isa Float32
+    @test pdf(LogUniform(1, 2), 1) isa Float64
+    @test quantile(LogUniform(1, 2), 1) isa Float64
+    @test quantile(LogUniform(1, 2), 1f0) isa Float32
+    @testset "$f" for f in [pdf, cdf, quantile, logpdf, logcdf]
+        @test @inferred(f(LogUniform(1,2), 1)) isa Float64
+        @test @inferred(f(LogUniform(1,2), 1.0)) isa Float64
+        @test @inferred(f(LogUniform(1.0,2), 1.0)) isa Float64
+        @test @inferred(f(LogUniform(1.0f0,2), 1)) isa Float32
+        @test @inferred(f(LogUniform(1.0f0,2), 1f0)) isa Float32
+        @test @inferred(f(LogUniform(1,2), 1f0)) isa Float32
+    end
+
+    d = LogUniform(1,10)
+    @test eltype(d) === Float64
+    @test 1 <= rand(rng, d) <= 10
+    @test rand(rng, d) isa eltype(d)
+    @test @inferred(quantile(d, 0))   ≈ 1
+    @test quantile(d, 0.5) ≈ sqrt(10) # geomean
+    @test quantile(d, 1)   ≈ 10
+    @test mode(d) ≈ 1
+    @test !insupport(d, 0)
+    @test @inferred(minimum(d)) === 1
+    @test @inferred(maximum(d)) === 10
+    @test partype(d) === Int
+    @test truncated(d, 2, 14) === LogUniform(2,10)
+
+    # numbers obtained by calling scipy.stats.loguniform
+    @test @inferred(std(d)        ) ≈ 2.49399867607628
+    @test @inferred(mean(d)       ) ≈ 3.908650337129266
+    @test @inferred(pdf(d, 1.0001)) ≈ 0.43425105679757203
+    @test @inferred(pdf(d, 5     )) ≈ 0.08685889638065035
+    @test @inferred(pdf(d, 9.9999)) ≈ 0.04342988248915007
+    @test @inferred(cdf(d, 1.0001)) ≈ 4.342727686266485e-05
+    @test @inferred(cdf(d, 5     )) ≈ 0.6989700043360187
+    @test @inferred(cdf(d, 9.9999)) ≈ 0.999995657033466
+    @test @inferred(median(d)     ) ≈ 3.1622776601683795
+    @test @inferred(logpdf(d, 5)  ) ≈ -2.443470357682056
+
+    for _ in 1:10
+        lo = rand(rng)
+        hi = lo + 10*rand(rng)
+        dist = LogUniform(lo,hi)
+        q = rand(rng)
+        @test cdf(dist, quantile(dist, q)) ≈ q
+
+        u = Uniform(log(lo), log(hi))
+        @test exp(quantile(u, q)) ≈ quantile(dist, q)
+        @test exp(median(u)) ≈ median(dist)
+        x = rand(rng, dist)
+        @test cdf(u, log(x)) ≈ cdf(dist, x)
+
+        @test @inferred(entropy(dist)) ≈  Distributions.expectation(dist, x->-logpdf(dist,x))
+    end
+
+    @test kldivergence(LogUniform(1,2), LogUniform(1,2)) ≈ 0 atol=100eps(Float64)
+    @test isfinite(kldivergence(LogUniform(1,2), LogUniform(1,10)))
+    @test kldivergence(LogUniform(1.1,10), LogUniform(1,2)) === Inf
+    @test kldivergence(LogUniform(0.1,10), LogUniform(1,2)) === Inf
+    @test kldivergence(LogUniform(0.1,1), LogUniform(1,2))  === Inf
+    @test @inferred(kldivergence(LogUniform(0.1f0,1), LogUniform(1,2)))  === Inf32
+
+    for _ in 1:10
+        aq = 10*rand(rng)
+        ap = aq + 10*rand(rng)
+        bp = ap + 10*rand(rng)
+        bq = bp + 10*rand(rng)
+        p = LogUniform(ap, bp)
+        q = LogUniform(aq, bq)
+        @test @inferred(kldivergence(p, q)) ≈
+            kldivergence(Uniform(log(ap), log(bp)), Uniform(log(aq), log(bq)))
+    end
+end
+
+
+end#module

test/runtests.jl

@@ -11,6 +11,7 @@ import JSON
 import ForwardDiff
 
 const tests = [
+    "loguniform",
     "arcsine",
     "dirac",
     "truncate",