Add the one-parameter Lindley distribution

.github/workflows/CI.yml

@@ -24,7 +24,7 @@ jobs:
       fail-fast: false
       matrix:
         version:
-          - '1.3'
+          - '1.6'
           - '1'
           - 'nightly'
         os:

Project.toml

@@ -7,6 +7,7 @@ version = "0.25.81"
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
 DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+LambertW = "984bce1d-4616-540c-a9ee-88d1112d94c9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PDMats = "90014a1f-27ba-587c-ab20-58faa44d9150"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -28,7 +29,7 @@ QuadGK = "2"
 SpecialFunctions = "1.2, 2"
 StatsBase = "0.32, 0.33"
 StatsFuns = "0.9.15, 1"
-julia = "1.3"
+julia = "1.6"
 
 [extras]
 Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"

docs/src/univariate.md

@@ -298,6 +298,13 @@ Levy
 plotdensity((0, 20), Levy, (0, 1)) # hide
 ```
 
+```@docs
+Lindley
+```
+```@example plotdensity
+plotdensity((0, 20), Lindley, (1.5,)) # hide
+```
+
 ```@docs
 Logistic
 ```

src/Distributions.jl

@@ -23,7 +23,7 @@ import StatsBase: kurtosis, skewness, entropy, mode, modes,
 
 import PDMats: dim, PDMat, invquad
 
-using SpecialFunctions
+using SpecialFunctions, LambertW
 
 import ChainRulesCore
 
@@ -114,6 +114,7 @@ export
     KSOneSided,
     Laplace,
     Levy,
+    Lindley,
     LKJ,
     LKJCholesky,
     LocationScale,
@@ -345,7 +346,7 @@ Supported distributions:
     Frechet, FullNormal, FullNormalCanon, Gamma, GeneralizedPareto,
     GeneralizedExtremeValue, Geometric, Gumbel, Hypergeometric,
     InverseWishart, InverseGamma, InverseGaussian, IsoNormal,
-    IsoNormalCanon, Kolmogorov, KSDist, KSOneSided, Laplace, Levy, LKJ, LKJCholesky,
+    IsoNormalCanon, Kolmogorov, KSDist, KSOneSided, Laplace, Levy, Lindley, LKJ, LKJCholesky,
     Logistic, LogNormal, MatrixBeta, MatrixFDist, MatrixNormal,
     MatrixTDist, MixtureModel, Multinomial,
     MultivariateNormal, MvLogNormal, MvNormal, MvNormalCanon,

src/univariate/continuous/lindley.jl

@@ -0,0 +1,184 @@
+"""
+    Lindley(θ)
+
+The one-parameter *Lindley distribution* with shape `θ > 0` has probability density
+function
+
+```math
+f(x; \\theta) = \\frac{\\theta^2}{1 + \\theta} (1 + x) e^{-\\theta x}, \\quad x > 0
+```
+
+It was first described by Lindley[^1] and was studied in greater detail by Ghitany
+et al.[^2]
+Note that `Lindley(θ)` is a mixture of an `Exponential(θ)` and a `Gamma(2, θ)` with
+respective mixing weights `p = θ/(1 + θ)` and `1 - p`.
+
+[^1]: Lindley, D. V. (1958). Fiducial Distributions and Bayes' Theorem. Journal of the
+      Royal Statistical Society: Series B (Methodological), 20(1), 102–107.
+[^2]: Ghitany, M. E., Atieh, B., & Nadarajah, S. (2008). Lindley distribution and its
+      application. Mathematics and Computers in Simulation, 78(4), 493–506.
+"""
+struct Lindley{T<:Real} <: ContinuousUnivariateDistribution
+    θ::T
+
+    Lindley{T}(θ::T) where {T} = new{T}(θ)
+end
+
+function Lindley(θ::Real; check_args::Bool=true)
+    @check_args Lindley (θ, θ > zero(θ))
+    return Lindley{typeof(θ)}(θ)
+end
+
+Lindley(θ::Integer; check_args::Bool=true) = Lindley(float(θ); check_args=check_args)
+
+Lindley() = Lindley{Float64}(1.0)
+
+Base.convert(::Type{Lindley{T}}, d::Lindley) where {T} = Lindley{T}(T(shape(d)))
+
+@distr_support Lindley 0.0 Inf
+
+### Parameters
+
+shape(d::Lindley) = d.θ
+params(d::Lindley) = (shape(d),)
+partype(::Lindley{T}) where {T} = T
+
+### Statistics
+
+mean(d::Lindley) = (2 + d.θ) / d.θ / (1 + d.θ)
+
+var(d::Lindley) = 2 / d.θ^2 - 1 / (1 + d.θ)^2
+
+skewness(d::Lindley) = 2 * evalpoly(d.θ, (2, 6, 6, 1)) / evalpoly(d.θ, (2, 4, 1))^(3//2)
+
+kurtosis(d::Lindley) = 3 * evalpoly(d.θ, (8, 32, 44, 24, 3)) / evalpoly(d.θ, (2, 4, 1))^2 - 3
+
+mode(d::Lindley) = d.θ < 1 ? (1 - d.θ) / d.θ : zero(d.θ)
+
+# Derived with Mathematica:
+#     KLDivergence := ResourceFunction["KullbackLeiblerDivergence"]
+#     KLDivergence[LindleyDistribution[θp], LindleyDistribution[θq]]
+function kldivergence(p::Lindley, q::Lindley)
+    θp = shape(p)
+    θq = shape(q)
+    a = (θp + 2) * (θp - θq) / θp / (1 + θp)
+    b = 2 * log(θp) + log1p(θq) - 2 * log(θq) - log1p(θp)
+    return b - a
+end
+
+# Derived with Mathematica based on https://mathematica.stackexchange.com/a/275765:
+#     ShannonEntropy[dist_?DistributionParameterQ] :=
+#         Expectation[-LogLikelihood[dist, {x}], Distributed[x, dist]]
+#     Simplify[ShannonEntropy[LindleyDistribution[θ]]]
+function entropy(d::Lindley)
+    θ = shape(d)
+    return 1 + exp(θ) * expinti(-θ) / (1 + θ) - 2 * log(θ) + log1p(θ)
+end
+
+### Evaluation
+
+_lindley_mgf(θ, t) = θ^2 * (1 + θ - t) / (1 + θ) / (θ - t)^2
+
+mgf(d::Lindley, t::Real) = _lindley_mgf(shape(d), t)
+
+cf(d::Lindley, t::Real) = _lindley_mgf(shape(d), t * im)
+
+cgf(d::Lindley, t::Real) = log1p(-t / (1 + d.θ)) - 2 * log1p(-t / d.θ)
+
+_zero(d::Lindley, y::Real) = zero(shape(d)) * zero(y)
+_oftype(d::Lindley, y::Real, x) = oftype(_zero(d, y), x)
+
+function pdf(d::Lindley, y::Real)
+    θ = shape(d)
+    if isnan(y)
+        return _oftype(d, y, NaN)
+    elseif isfinite(y) && y > 0
+        return θ^2 / (1 + θ) * (1 + y) * exp(-θ * y)
+    else
+        return _zero(d, y)
+    end
+end
+
+function logpdf(d::Lindley, y::Real)
+    θ = shape(d)
+    if isnan(y)
+        return _oftype(d, y, NaN)
+    elseif isfinite(y) && y > 0
+        return 2 * log(θ) - log1p(θ) + log1p(y) - θ * y
+    else
+        return _oftype(d, y, -Inf)
+    end
+end
+
+function gradlogpdf(d::Lindley, y::Real)
+    if isnan(y)
+        return _oftype(d, y, NaN)
+    elseif isfinite(y) && y > 0
+        return inv(1 + y) - shape(d)
+    else
+        return _zero(d, y)
+    end
+end
+
+function ccdf(d::Lindley, y::Real)
+    θ = shape(d)
+    θy = θ * y
+    if isnan(y)
+        return _oftype(d, y, NaN)
+    elseif y > 0
+        if isfinite(y)
+            return (1 + θy / (1 + θ)) * exp(-θy)
+        else
+            return _zero(d, y)
+        end
+    else
+        return _oftype(d, y, 1)
+    end
+end
+
+function logccdf(d::Lindley, y::Real)
+    θ = shape(d)
+    if isnan(y)
+        return _oftype(d, y, NaN)
+    elseif y > 0
+        if isfinite(y)
+            return log1p(θ * (1 + y)) - log1p(θ) - θ * y
+        else
+            return _oftype(d, y, -Inf)
+        end
+    else
+        return _zero(d, y)
+    end
+end
+
+cdf(d::Lindley, y::Real) = 1 - ccdf(d, y)
+
+logcdf(d::Lindley, y::Real) = log1p(-ccdf(d, y))
+
+# Jodrá, P. (2010). Computer generation of random variables with Lindley or
+# Poisson–Lindley distribution via the Lambert W function. Mathematics and Computers
+# in Simulation, 81(4), 851–859.
+function quantile(d::Lindley, q::Real)
+    θ = shape(d)
+    return -1 - inv(θ) - lambertw((1 + θ) * (q - 1) / exp(1 + θ), -1) / θ
+end
+
+### Sampling
+
+# Ghitany, M. E., Atieh, B., & Nadarajah, S. (2008). Lindley distribution and its
+# application. Mathematics and Computers in Simulation, 78(4), 493–506.
+function rand(rng::AbstractRNG, d::Lindley{T}) where {T}
+    θ = shape(d)
+    λ = inv(θ)
+    u = rand(rng, T)
+    p = θ / (1 + θ)
+    return T(rand(rng, u <= p ? Exponential{T}(λ) : Gamma{T}(2, λ)))
+end
+
+### Fitting
+
+# Ghitany et al. (2008)
+function fit_mle(::Type{<:Lindley}, x::AbstractArray{<:Real})
+    x̄ = mean(x)
+    return Lindley((1 - x̄ + sqrt((x̄ - 1)^2 + 8x̄)) / 2x̄)
+end

src/univariates.jl

@@ -692,6 +692,7 @@ const continuous_distributions = [
     "ksonesided",
     "laplace",
     "levy",
+    "lindley",
     "logistic",
     "noncentralbeta",
     "noncentralchisq",

test/ref/continuous/lindley.R

@@ -0,0 +1,12 @@
+library("LindleyR")
+
+Lindley <- R6Class("Lindley",
+                   inherit=ContinuousDistribution,
+                   public=list(names=c("theta"),
+                               theta=NA,
+                               initialize=function(theta=1) { self$theta <- theta },
+                               supp=function() { c(0, Inf) },
+                               properties=function() { list() },
+                               pdf=function(x, log=FALSE) { dlindley(x, self$theta, log=log) },
+                               cdf=function(x) { plindley(x, self$theta) },
+                               quan=function(x) { qlindley(x, self$theta) }))

test/ref/continuous_test.lst

@@ -101,6 +101,11 @@ Levy(2)
 Levy(2, 8)
 Levy(3.0, 3)
 
+Lindley()
+Lindley(0.5)
+Lindley(1.5)
+Lindley(3.0)
+
 Logistic()
 Logistic(2.0)
 Logistic(0.0, 1.0)

test/ref/continuous_test.ref.json

@@ -2876,6 +2876,110 @@
     { "q": 0.90, "x": 2.19722457733622 }
   ]
 },
+{
+  "expr": "Lindley()",
+  "dtype": "Lindley",
+  "minimum": 0,
+  "maximum": "inf",
+  "properties": {
+  },
+  "points": [
+    { "x": 0.201229322235454, "pdf": 0.491137556224262, "logpdf": -0.711031035180873, "cdf": 0.1 },
+    { "x": 0.409355765160304, "pdf": 0.467961032750689, "logpdf": -0.759370249884304, "cdf": 0.2 },
+    { "x": 0.630825401024496, "pdf": 0.433923809718651, "logpdf": -0.834886313936738, "cdf": 0.3 },
+    { "x": 0.873054028699824, "pdf": 0.391162994497697, "logpdf": -0.938630940138138, "cdf": 0.4 },
+    { "x": 1.14619322062058, "pdf": 0.341077783550314, "logpdf": -1.07564472345732, "cdf": 0.5 },
+    { "x": 1.4662034354858, "pdf": 0.284599964357274, "logpdf": -1.25667071856436, "cdf": 0.6 },
+    { "x": 1.86201458406329, "pdf": 0.222320334770884, "logpdf": -1.5036359877426, "cdf": 0.7 },
+    { "x": 2.39727598728087, "pdf": 0.154517296485711, "logpdf": -1.8674492375466, "cdf": 0.8 },
+    { "x": 3.27181206035629, "pdf": 0.0810311902520209, "logpdf": -2.51292113358837, "cdf": 0.9 }
+  ],
+  "quans": [
+    { "q": 0.10, "x": 0.201229322235454 },
+    { "q": 0.25, "x": 0.517999713886834 },
+    { "q": 0.50, "x": 1.14619322062058 },
+    { "q": 0.75, "x": 2.10546657787674 },
+    { "q": 0.90, "x": 3.27181206035629 }
+  ]
+},
+{
+  "expr": "Lindley(0.5)",
+  "dtype": "Lindley",
+  "minimum": 0,
+  "maximum": "inf",
+  "properties": {
+  },
+  "points": [
+    { "x": 0.544019554971776, "pdf": 0.196051062631011, "logpdf": -1.6293801300544, "cdf": 0.0999999999999999 },
+    { "x": 1.04307224421496, "pdf": 0.202130669036478, "logpdf": -1.59884091429709, "cdf": 0.2 },
+    { "x": 1.5435367816863, "pdf": 0.195934998738978, "logpdf": -1.62997231384287, "cdf": 0.3 },
+    { "x": 2.07183011333304, "pdf": 0.181699507555922, "logpdf": -1.70540101378774, "cdf": 0.4 },
+    { "x": 2.65368480453801, "pdf": 0.161562102012007, "logpdf": -1.82286567765105, "cdf": 0.5 },
+    { "x": 3.32408881228643, "pdf": 0.136749781372001, "logpdf": -1.98960243642355, "cdf": 0.6 },
+    { "x": 4.14298158855248, "pdf": 0.10800073172794, "logpdf": -2.22561727662217, "cdf": 0.7 },
+    { "x": 5.23954037209116, "pdf": 0.0757268013786986, "logpdf": -2.58062313392393, "cdf": 0.8 },
+    { "x": 7.01639138849524, "pdf": 0.0400163645647015, "logpdf": -3.21846679441503, "cdf": 0.9 }
+  ],
+  "quans": [
+    { "q": 0.10, "x": 0.544019554971776 },
+    { "q": 0.25, "x": 1.29133558854642 },
+    { "q": 0.50, "x": 2.65368480453801 },
+    { "q": 0.75, "x": 4.64292488017037 },
+    { "q": 0.90, "x": 7.01639138849524 }
+  ]
+},
+{
+  "expr": "Lindley(1.5)",
+  "dtype": "Lindley",
+  "minimum": 0,
+  "maximum": "inf",
+  "properties": {
+  },
+  "points": [
+    { "x": 0.114556971437579, "pdf": 0.844729364271255, "logpdf": -0.168738981894155, "cdf": 0.0999999999999998 },
+    { "x": 0.237615294552685, "pdf": 0.77989414577674, "logpdf": -0.248597079050649, "cdf": 0.2 },
+    { "x": 0.372143951935997, "pdf": 0.706662569042456, "logpdf": -0.347201998525827, "cdf": 0.3 },
+    { "x": 0.522280435214869, "pdf": 0.625895614615692, "logpdf": -0.468571671601523, "cdf": 0.4 },
+    { "x": 0.694246765710091, "pdf": 0.538217563107918, "logpdf": -0.619492408171021, "cdf": 0.5 },
+    { "x": 0.89826450068756, "pdf": 0.444050395936119, "logpdf": -0.811817218630417, "cdf": 0.6 },
+    { "x": 1.15323211442359, "pdf": 0.343613202710769, "logpdf": -1.06823866495725, "cdf": 0.7 },
+    { "x": 1.50109080206919, "pdf": 0.236863853380854, "logpdf": -1.44026976122003, "cdf": 0.8 },
+    { "x": 2.07401891659104, "pdf": 0.12326693255173, "logpdf": -2.09340309170877, "cdf": 0.9 }
+  ],
+  "quans": [
+    { "q": 0.10, "x": 0.114556971437579 },
+    { "q": 0.25, "x": 0.30322247911811 },
+    { "q": 0.50, "x": 0.694246765710091 },
+    { "q": 0.75, "x": 1.31109339979007 },
+    { "q": 0.90, "x": 2.07401891659104 }
+  ]
+},
+{
+  "expr": "Lindley(3.0)",
+  "dtype": "Lindley",
+  "minimum": 0,
+  "maximum": "inf",
+  "properties": {
+  },
+  "points": [
+    { "x": 0.0465620413518296, "pdf": 2.04777663835177, "logpdf": 0.716754637924627, "cdf": 0.1 },
+    { "x": 0.0980294998228117, "pdf": 1.84109209665868, "logpdf": 0.610358926343917, "cdf": 0.2 },
+    { "x": 0.155708594360887, "pdf": 1.62989906665426, "logpdf": 0.488518090603404, "cdf": 0.3 },
+    { "x": 0.221505053392774, "pdf": 1.41410780366481, "logpdf": 0.34649880478381, "cdf": 0.4 },
+    { "x": 0.298361440425077, "pdf": 1.19357014066527, "logpdf": 0.176948933955783, "cdf": 0.5 },
+    { "x": 0.391185596605387, "pdf": 0.968051255886061, "logpdf": -0.0324702428118596, "cdf": 0.6 },
+    { "x": 0.509131699567543, "pdf": 0.737174657514306, "logpdf": -0.304930430453559, "cdf": 0.7 },
+    { "x": 0.672626102659318, "pdf": 0.500297086565447, "logpdf": -0.692553183880015, "cdf": 0.8 },
+    { "x": 0.94630688386014, "pdf": 0.256133428755212, "logpdf": -1.36205676420823, "cdf": 0.9 }
+  ],
+  "quans": [
+    { "q": 0.10, "x": 0.0465620413518296 },
+    { "q": 0.25, "x": 0.125991234674835 },
+    { "q": 0.50, "x": 0.298361440425077 },
+    { "q": 0.75, "x": 0.583010462004234 },
+    { "q": 0.90, "x": 0.94630688386014 }
+  ]
+},
 {
   "expr": "Logistic(2.0)",
   "dtype": "Logistic",