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,142 @@
+"""
+    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{Logistic}(\\theta, \\phi)`` then ``exp(X) \\sim \\operatorname{LogLogistic}(\\theta, \\phi)``. The probability density function is 
+
+```math
+f(x; \\theta, \\phi) = \\frac{(\\phi / \\theta)x/\\theta()^(\\phi - 1)}{(1 + (x/\\theta)^\\phi)^2}, \\theta > 0, \\phi > 0
+```
+
+```julia
+LogLogistic()            # Log-logistic distribution with unit scale and unit shape
+LogLogistic(θ)           # Log-logistic distribution with 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(θ::T) where {T<:Real} = LogLogistic(θ, 1.0)
+LogLogistic() = LogLogistic(1.0, 1.0, check_args=false) 
+
+@distr_support LogLogistic 0.0 Inf
+
+#### Coversions
+convert(::Type{LogLogistic{T}}, θ::S, ϕ::S) where {T <: Real, S <: Real} = LogLogistic(T(θ), T(ϕ))
+convert(::Type{LogLogistic{T}}, d::LogLogistic{S}) where {T <: Real, S <: Real} = LogLogistic(T(d.θ), T(d.ϕ), check_args=false)
+
+#### Parameters 
+
+params(d::LogLogistic) = (d.θ, d.ϕ)
+partype(::LogLogistic{T}) where {T} = T
+
+#### Statistics 
+
+median(d::LogLogistic) = d.θ
+function mean(d::LogLogistic)
+	if d.ϕ ≤ 1
+        error("mean is defined only when ϕ > 1") 	
+	end
+	return d.θ*π/d.ϕ/sin(π/d.ϕ)
+end
+
+function mode(d::LogLogistic)
+	if d.ϕ ≤ 1
+		error("mode is defined only when ϕ > 1")
+	end 
+	return d.θ*((d.ϕ-1)/(d.ϕ+1))^(1/d.ϕ)
+end
+
+function var(d::LogLogistic)
+	if d.ϕ ≤ 2
+		erros("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, x::Real)
+    if x ≤ zero(0)
+        z = zero(x)
+    else
+        # use built-in impletation to evaluate the density 
+        # of loglogistic at x 
+        # Y = log(X)
+        # Y ~ logistic(log(θ), 1/ϕ)
+        z = pdf(Logistic(log(d.θ), 1/d.ϕ), log(x)) / x
+    end    
+    return z
+end
+
+function logpdf(d::LogLogistic, x::Real)
+    if x ≤ zero(0)
+        z = log(zero(x))
+    else
+        z = logpdf(Logistic(log(d.θ), 1/d.ϕ), log(x)) + log(x)
+    end
+    return z 
+end
+
+function cdf(d::LogLogistic, x::Real)
+    if x <= 0
+        return 0.0
+    end 
+    z = cdf(Logistic(log(d.θ), 1/d.ϕ), log(x)) 
+    return z
+end
+
+function logcdf(d::LogLogistic, x::Real)
+    if x <= 0
+        -Inf
+    end
+    z = logcdf(Logistic(log(d.θ), 1/d.ϕ), log(x)) 
+	return z
+end
+
+function ccdf(d::LogLogistic, x::Real)
+    if x <= 0
+        return 1
+    end 
+    z = ccdf(Logistic(log(d.θ), 1/d.ϕ), log(x)) 
+	return z 
+end
+
+function logccdf(d::LogLogistic, x::Real)
+    if x <= 0
+        return 0.0
+    end
+    z = logccdf(Logistic(log(d.θ), 1/d.ϕ), log(x)) 
+	return z 
+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,25 @@
+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), digits=6) == 0.125
+    @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), 3), digits=6) == 0.6
+    @test round(cdf(LogLogistic(2,2), 4), digits=6) == 0.8
+
+    @test round(ccdf(LogLogistic(2), 3), digits=6) == 0.4
+    @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