Add log-logistic distribution

Project.toml

@@ -1,7 +1,7 @@
 name = "Distributions"
 uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
 authors = ["JuliaStats"]
-version = "0.25.122"
+version = "0.25.123"
 
 [deps]
 AliasTables = "66dad0bd-aa9a-41b7-9441-69ab47430ed8"
@@ -43,6 +43,7 @@ FiniteDifferences = "0.12"
 ForwardDiff = "0.10, 1"
 JSON = "0.21"
 LinearAlgebra = "<0.0.1, 1"
+Optim = "1.13"
 OffsetArrays = "1"
 PDMats = "0.11.35"
 Printf = "<0.0.1, 1"
@@ -69,11 +70,12 @@ Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
 FiniteDifferences = "26cc04aa-876d-5657-8c51-4c34ba976000"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
+Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "StableRNGs", "Calculus", "ChainRulesCore", "ChainRulesTestUtils", "DensityInterface", "Distributed", "FiniteDifferences", "ForwardDiff", "JSON", "SparseArrays", "StaticArrays", "Test", "OffsetArrays"]
+test = ["Aqua", "StableRNGs", "Calculus", "ChainRulesCore", "ChainRulesTestUtils", "DensityInterface", "Distributed", "FiniteDifferences", "ForwardDiff", "JSON", "Optim", "SparseArrays", "StaticArrays", "Test", "OffsetArrays"]

docs/Project.toml

@@ -1,7 +1,11 @@
 [deps]
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
 
 [compat]
 Documenter = "1"
 GR = "0.72.1, 0.73"
+
+[sources.Distributions]
+path = ".."

docs/src/univariate.md

@@ -326,6 +326,13 @@ Logistic
 plotdensity((-4, 8), Logistic, (2, 1)) # hide
 ```
 
+```@docs
+LogLogistic
+```
+```@example plotdensity
+plotdensity((0, 2), LogLogistic, (2, 1)) # hide
+```
+
 ```@docs
 LogitNormal
 ```

src/Distributions.jl

@@ -123,6 +123,7 @@ export
     LKJCholesky,
     LocationScale,
     Logistic,
+    LogLogistic,
     LogNormal,
     LogUniform,
     MvLogitNormal,
@@ -354,7 +355,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,120 @@
+"""
+    LogLogistic(α, β)
+
+The *log-logistic distribution* with scale ``\\alpha`` and shape ``\\beta`` 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{(\\beta / \\alpha){(x/\\alpha)}^{\\beta - 1}}{{(1 + {(x/\\alpha)}^\\beta)}^2}, \\qquad \\alpha > 0, \\beta > 0.
+```
+
+```julia
+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::Bool=true) where {T <: Real}
+    check_args && @check_args(LogLogistic, α > zero(α) && β > zero(β))
+    return LogLogistic{T}(α, β)
+end
+
+LogLogistic(α::Real, β::Real; check_args::Bool = true) = LogLogistic(promote(α, β)...; check_args)
+
+@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)
+    (; α, β) = d
+	if !(β > 1)
+        throw(ArgumentError("the mean of a log-logistic distribution is defined only when its shape β > 1"))
+	end
+    return α/sinc(inv(β))
+end
+
+function mode(d::LogLogistic)
+    (; α, β) = d
+    return α*(max(β - 1, 0) / (β + 1))^inv(β)
+end
+
+function var(d::LogLogistic)
+    (; α, β) = d
+	if !(β > 2)
+        throw(ArgumentError("the variance of a log-logistic distribution is defined only when its shape β > 2"))
+	end
+    invβ = inv(β)
+	return α^2 * (inv(sinc(2 * invβ)) - inv(sinc(invβ))^2)
+end
+
+entropy(d::LogLogistic) = log(d.α / d.β) + 2
+
+#### Evaluation
+
+function pdf(d::LogLogistic, x::Real)
+    (; α, β) = d
+    insupport = x > 0
+    _x = insupport ? x : zero(x)
+    xoαβ = (_x / α)^β
+    res = (β / _x) / ((1 + xoαβ) * (1 + inv(xoαβ)))
+    return insupport ? res : zero(res)
+end
+function logpdf(d::LogLogistic, x::Real)
+    (; α, β) = d
+    insupport = x > 0
+    _x = insupport ? x : zero(x)
+    βlogxoα = β * log(_x / α)
+    res = log(β / _x) - (log1pexp(βlogxoα) + log1pexp(-βlogxoα))
+    return insupport ? res : oftype(res, -Inf)
+end
+
+cdf(d::LogLogistic, x::Real) = inv(1 + (max(x, 0) / d.α)^(-d.β))
+ccdf(d::LogLogistic, x::Real) = inv(1 + (max(x, 0) / d.α)^d.β)
+
+logcdf(d::LogLogistic, x::Real) = -log1pexp(-d.β * log(max(x, 0) / d.α))
+logccdf(d::LogLogistic, x::Real) = -log1pexp(d.β * log(max(x, 0) / d.α))
+
+quantile(d::LogLogistic, p::Real) = d.α * (p / (1 - p))^inv(d.β)
+cquantile(d::LogLogistic, p::Real) = d.α * ((1 - p) / p)^inv(d.β)
+
+invlogcdf(d::LogLogistic, lp::Real) = d.α * expm1(-lp)^(-inv(d.β))
+invlogccdf(d::LogLogistic, lp::Real) = d.α * expm1(-lp)^inv(d.β)
+
+#### Sampling
+
+function rand(rng::AbstractRNG, d::LogLogistic)
+    u = rand(rng)
+    return d.α * (u / (1 - u))^(inv(d.β))
+end
+function rand!(rng::AbstractRNG, d::LogLogistic, A::AbstractArray{<:Real})
+    rand!(rng, A)
+    let α = d.α, invβ = inv(d.β)
+        map!(A, A) do u
+            return α * (u / (1 - u))^invβ
+        end
+    end
+    return A
+end

src/univariates.jl

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

test/ref/continuous/loglogistic.R

@@ -0,0 +1,17 @@
+LogLogistic <- R6Class("LogLogistic",
+    inherit = ContinuousDistribution,
+    public = list(
+        names = c("alpha", "beta"),
+        alpha = NA,
+        beta = NA,
+        initialize = function(a, b) {
+            self$alpha <- a
+            self$beta <- b
+        },
+        supp = function() { c(0, Inf) },
+        properties = function() { },
+        pdf = function(x, log=FALSE) { VGAM::dfisk(x, scale=self$alpha, shape=self$beta, log=log) },
+        cdf = function(x) { VGAM::pfisk(x, scale=self$alpha, shape=self$beta) },
+        quan = function(v) { VGAM::qfisk(v, scale=self$alpha, shape=self$beta) }
+    )
+)

test/ref/continuous_test.lst

@@ -126,6 +126,14 @@ Logistic(5.0, 1.0)
 Logistic(2.0, 1.5)
 Logistic(5.0, 1.5)
 
+LogLogistic(0.5, 0.5)
+LogLogistic(0.5, 1.0)
+LogLogistic(1.0, 0.5)
+LogLogistic(1.0, 1.0)
+LogLogistic(1.0, 2.0)
+LogLogistic(2.0, 1.0)
+LogLogistic(2.0, 2.0)
+
 LogNormal()
 LogNormal(1.0)
 LogNormal(0.0, 2.0)