add loglogistic distribution

docs/Project.toml

@@ -1,4 +1,5 @@
 [deps]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
 

docs/src/univariate.md

@@ -320,10 +320,17 @@ plotdensity((0, 20), Lindley, (1.5,)) # hide
 ```
 
 ```@docs
-Logistic
+logistic
 ```
 ```@example plotdensity
-plotdensity((-4, 8), Logistic, (2, 1)) # hide
+plotdensity((-4, 8), logistic, (2, 1)) # hide
+```
+
+```@docs
+loglogistic
+```
+```@example plotdensity
+plotdensity((0, 2), logistic, (1, 1)) # hide
 ```
 
 ```@docs

src/Distributions.jl

@@ -120,6 +120,7 @@ export
     LKJCholesky,
     LocationScale,
     Logistic,
+    LogLogistic,
     LogNormal,
     LogUniform,
     LogitNormal,
@@ -356,7 +357,7 @@ Supported distributions:
     InverseWishart, InverseGamma, InverseGaussian, IsoNormal,
     IsoNormalCanon, JohnsonSU, Kolmogorov, KSDist, KSOneSided, Kumaraswamy,
     Laplace, Levy, Lindley, LKJ, LKJCholesky,
-    Logistic, LogNormal, MatrixBeta, MatrixFDist, MatrixNormal,
+    Logistic, LogLogistic, LogNormal, MatrixBeta, MatrixFDist, MatrixNormal,
     MatrixTDist, MixtureModel, Multinomial,
     MultivariateNormal, MvLogNormal, MvNormal, MvNormalCanon,
     MvNormalKnownCov, MvTDist, NegativeBinomial, NoncentralBeta, NoncentralChisq,

src/univariate/continuous/loglogistic.jl

@@ -0,0 +1,111 @@
+"""
+    LogLogistic(α, β)
+
+The *log logistic distribution* with scale `α` and shape `β` is the distribution of a random variable whose logarithm has a [`Logistic`](@ref) distribution. 
+If ``X \\sim \\operatorname{LogLogistic}(\\alpha, \\beta)`` then ``log(X) \\sim \\operatorname{Logistic}(log(\\alpha), 1/\\beta)``. The probability density function is 
+
+```math
+f(x; \\alpha, \\beta) = \\frac{(\\alpha / \\beta)x/\\beta()^(\\alpha - 1)}{(1 + (x/\\beta)^\\alpha)^2}, \\beta > 0, \\alpha > 0
+```
+
+```julia
+LogLogistic()            # Log-logistic distribution with unit scale and unit shape
+LogLogistic(α,β)         # Log-logistic distribution with scale α and shape β
+
+params(d)                # Get the parameters, i.e. (α, β)
+scale(d)                 # Get the scale parameter, i.e. α
+shape(d)                 # Get the shape parameter, i.e. β
+```
+
+External links
+
+* [Log logistic distribution on Wikipedia](https://en.wikipedia.org/wiki/Log-logistic_distribution)
+"""
+
+struct LogLogistic{T<:Real} <: ContinuousUnivariateDistribution
+    α::T
+    β::T
+    LogLogistic{T}(α::T,β::T) where {T} = new{T}(α,β)
+end
+
+function LogLogistic(α::T, β::T; check_args=true) where {T <: Real}
+    check_args && @check_args(LogLogistic, α > zero(α) && β > zero(β))
+    return LogLogistic{T}(α, β)
+end
+
+LogLogistic(α::Real, β::Real) = LogLogistic(promote(α,β)...)
+LogLogistic(α::Integer, β::Integer) = LogLogistic(float(α), float(β))
+LogLogistic() = LogLogistic(1.0, 1.0, check_args=false) 
+
+@distr_support LogLogistic 0.0 Inf
+
+#### Coversions
+convert(::Type{LogLogistic{T}}, d::LogLogistic{T}) where {T<:Real} = d
+convert(::Type{LogLogistic{T}}, d::LogLogistic) where {T<:Real} = LogLogistic{T}(T(d.α), T(d.β))
+#### Parameters 
+
+params(d::LogLogistic) = (d.α, d.β)
+partype(::LogLogistic{T}) where {T} = T
+
+#### Statistics 
+
+median(d::LogLogistic) = d.α
+function mean(d::LogLogistic{T}) where T<:Real
+	if d.β ≤ 1
+        ArgumentError("mean is defined only when β > 1") 	
+	end
+    return d.α/sinc(1/d.β)
+end
+
+function mode(d::LogLogistic{T}) where T<:Real
+	if d.β ≤ 1
+		ArgumentError("mode is defined only when β > 1")
+	end 
+	return d.α*((d.β-1)/(d.β+1))^(1/d.β)
+end
+
+function var(d::LogLogistic{T}) where T<:Real
+	if d.β ≤ 2
+		ArgumentError("var is defined only when β > 2") 
+	end
+    b = π/d.β
+	return d.α^2 * (2*b/sin(2*b)-b^2/(sin(b))^2)
+end
+
+
+#### Evaluation
+function pdf(d::LogLogistic{T}, x::Real) where T<:Real
+    # use built-in impletation to evaluate the density 
+    # of loglogistic at x 
+    # Y = log(X)
+    # Y ~ logistic(log(θ), 1/ϕ)
+    x >= 0 ? pdf(Logistic(log(d.α), 1/d.β), log(x)) / x : zero(T)
+end
+
+function logpdf(d::LogLogistic{T}, x::Real) where T<:Real
+    x >= 0 ? logpdf(Logistic(log(d.α), 1/d.β), log(x)) + log(x) : -T(Inf)
+end
+
+function cdf(d::LogLogistic{T}, x::Real) where T<:Real
+    x >= 0 ? cdf(Logistic(log(d.α), 1/d.β), log(x)) : zero(T)
+end
+
+function logcdf(d::LogLogistic{T}, x::Real) where T<:Real
+    x >= 0 ? logcdf(Logistic(log(d.α), 1/d.β), log(x)) : -T(Inf)
+end
+
+function ccdf(d::LogLogistic{T}, x::Real) where T<:Real
+    x >= 0 ? ccdf(Logistic(log(d.α), 1/d.β), log(x)) : one(T)
+end
+
+function logccdf(d::LogLogistic{T}, x::Real) where T<:Real
+    x >= 0 ? logccdf(Logistic(log(d.α), 1/d.β), log(x)) : zero(T)
+end
+
+
+#### Sampling
+function rand(rng::AbstractRNG, d::LogLogistic)
+    u = rand(rng)
+    r = u / (1 - u)
+    return r^(1/d.β)*d.α
+end

src/univariates.jl

@@ -697,6 +697,7 @@ const continuous_distributions = [
     "levy",
     "lindley",
     "logistic",
+    "loglogistic",
     "noncentralbeta",
     "noncentralchisq",
     "noncentralf",

test/runtests.jl

@@ -11,6 +11,7 @@ import JSON
 import ForwardDiff
 
 const tests = [
+    "univariate/continuous/loglogistic",
     "univariate/continuous/loguniform",
     "univariate/continuous/arcsine",
     "univariate/discrete/dirac",
@@ -77,6 +78,7 @@ const tests = [
     "univariate/continuous/gumbel",
     "univariate/continuous/lindley",
     "univariate/continuous/logistic",
+    # "univariate/continuous/loglogistic",
     "univariate/continuous/johnsonsu",
     "univariate/continuous/noncentralchisq",
     "univariate/continuous/weibull",

test/univariate/continuous/loglogistic.jl

@@ -0,0 +1,22 @@
+using Distributions
+using Test
+
+@testset "LogLogistic" begin 
+
+    @test round(pdf(LogLogistic(), -1), digits=6) == 0.0
+    @test round(pdf(LogLogistic(), 1), digits=6) == 0.25
+    @test round(pdf(LogLogistic(2,2), 1), digits=6) == 0.32
+
+    @test round(cdf(LogLogistic(), -1), digits=6) == 0.0
+    @test round(cdf(LogLogistic(), 1), digits=6) == 0.5
+    @test round(cdf(LogLogistic(2,2), 4), digits=6) == 0.8
+
+    @test round(ccdf(LogLogistic(2,2), 4), digits=6) == 0.2
+
+    @test round(logpdf(LogLogistic(), -1), digits=6) == -Inf
+    @test round(logpdf(LogLogistic(), 1), digits=6) == -1.386294
+
+    @test round(logcdf(LogLogistic(2,2), 4), digits=6) == -0.223144
+    @test round(logccdf(LogLogistic(2,2), 4), digits=6) == -1.609438 
+
+end