Add OrderStatistic and JointOrderStatistics distributions

src/Distributions.jl

@@ -106,6 +106,7 @@ export
     InverseGaussian,
     IsoNormal,
     IsoNormalCanon,
+    JointOrderStatistics,
     Kolmogorov,
     KSDist,
     KSOneSided,
@@ -141,6 +142,7 @@ export
     Normal,
     NormalCanon,
     NormalInverseGaussian,
+    OrderStatistic,
     Pareto,
     PGeneralizedGaussian,
     SkewedExponentialPower,

src/multivariate/jointorderstatistics.jl

@@ -0,0 +1,132 @@
+# Implementation based on chapters 2-4 of
+# Arnold, Barry C., Narayanaswamy Balakrishnan, and Haikady Navada Nagaraja.
+# A first course in order statistics. Society for Industrial and Applied Mathematics, 2008.
+
+"""
+    JointOrderStatistics <: ContinuousMultivariateDistribution
+
+The joint distribution of a subset of order statistics from a sample from a continuous univariate distribution.
+
+    JointOrderStatistics(
+        dist::ContinuousUnivariateDistribution,
+        n::Int,
+        r::AbstractVector{Int}=1:n;
+        check_args::Bool=true,
+    )
+
+Construct the joint distribution of order statistics `r` from an IID sample of size `n` from `dist`.
+
+The ``i``th order statistic of a sample is the ``i``th element of the sorted sample.
+For example, the 1st order statistic is the sample minimum, while the ``n``th order
+statistic is the sample maximum.
+
+`r` must be a sorted vector of unique integers between 1 and `n`.
+
+For a single order statistic, use [`OrderStatistic`](@ref) instead.
+
+## Examples
+
+```julia
+JointOrderStatistics(Normal(), 10)           # Product(fill(Normal(), 10)) restricted to ordered vectors
+JointOrderStatistics(Cauchy(), 10, [1, 10])  # joint distribution of the extrema
+```
+"""
+struct JointOrderStatistics{D<:ContinuousUnivariateDistribution,R<:AbstractVector{Int}} <:
+       ContinuousMultivariateDistribution
+    dist::D
+    n::Int
+    r::R
+    function JointOrderStatistics(
+        dist::ContinuousUnivariateDistribution,
+        n::Int,
+        r::AbstractVector{Int}=1:n;
+        check_args::Bool=true,
+    )
+        @check_args(
+            JointOrderStatistics,
+            (n, n ≥ 1, "`n` must be a positive integer."),
+            (
+                r,
+                1 ≤ first(r) && last(r) ≤ n && issorted(r) && allunique(r),
+                "`r` must be a sorted vector of unique integers between 1 and `n`.",
+            ),
+        )
+        return new{typeof(dist),typeof(r)}(dist, n, r)
+    end
+end
+
+length(d::JointOrderStatistics) = length(d.r)
+function insupport(d::JointOrderStatistics, x::AbstractVector)
+    return length(d) == length(x) && issorted(x) && all(Base.Fix1(insupport, d.dist), x)
+end
+minimum(d::JointOrderStatistics) = fill(minimum(d.dist), length(d))
+maximum(d::JointOrderStatistics) = fill(maximum(d.dist), length(d))
+
+params(d::JointOrderStatistics) = tuple(params(d.dist)..., d.n, d.r)
+partype(d::JointOrderStatistics) = partype(d.dist)
+Base.eltype(::Type{<:JointOrderStatistics{D}}) where {D} = Base.eltype(D)
+
+function logpdf(d::JointOrderStatistics, x::AbstractVector{<:Real})
+    n = d.n
+    r = d.r
+    lp = sum(Base.Fix1(logpdf, d.dist), x)
+    T = eltype(lp)
+    lp += loggamma(T(n + 1))
+    length(r) == n && return lp
+    i = 0
+    xᵢ = oftype(float(first(x)), -Inf)
+    for (j, xⱼ) in zip(r, x)
+        lp += _marginalize_range(d.dist, n, i, j, xᵢ, xⱼ, T)
+        i = j
+        xᵢ = xⱼ
+    end
+    j = n + 1
+    xⱼ = oftype(xᵢ, Inf)
+    lp += _marginalize_range(d.dist, n, i, j, xᵢ, xⱼ, T)
+    return lp
+end
+
+# given ∏ₖf(xₖ), marginalize all xₖ for i < k < j
+function _marginalize_range(dist, n, i, j, xᵢ, xⱼ, T)
+    k = j - i - 1
+    k == 0 && return zero(T)
+    lpdiff = if i == 0
+        T(logcdf(dist, xⱼ))
+    elseif j == n + 1
+        T(logccdf(dist, xᵢ))
+    else
+        T(logdiffcdf(dist, xⱼ, xᵢ))
+    end
+    return k * lpdiff - loggamma(T(k + 1))
+end
+
+function _rand!(rng::AbstractRNG, d::JointOrderStatistics, x::AbstractVector{<:Real})
+    n = d.n
+    if n == length(d.r)  # r == 1:n
+        # direct method, slower than inversion method for large `n` and distributions with
+        # fast quantile function or that use inversion sampling
+        rand!(rng, d.dist, x)
+        sort!(x)
+    else
+        # use exponential generation method with inversion, where for gaps in the ranks, we
+        # use the fact that the sum Y of k IID variables xₘ ~ Exp(1) is Y ~ Gamma(k, 1).
+        # Lurie, D., and H. O. Hartley. "Machine-generation of order statistics for Monte
+        # Carlo computations." The American Statistician 26.1 (1972): 26-27.
+        # this is slow if length(d.r) is close to n and quantile for d.dist is expensive,
+        # but this branch is probably taken when length(d.r) is small or much smaller than n.
+        T = typeof(one(eltype(x)))
+        s = zero(one(T))
+        i = 0
+        for (m, j) in zip(eachindex(x), d.r)
+            k = j - i
+            s += k > 1 ? rand(rng, Gamma(k, one(T); check_args=false)) : randexp(rng, T)
+            i = j
+            x[m] = s
+        end
+        j = n + 1
+        k = j - i
+        s += k > 1 ? rand(rng, Gamma(k, one(T); check_args=false)) : randexp(rng, T)
+        x .= quantile.(d.dist, x ./ s)
+    end
+    return x
+end

src/multivariates.jl

@@ -112,6 +112,7 @@ end
 for fname in ["dirichlet.jl",
               "multinomial.jl",
               "dirichletmultinomial.jl",
+              "jointorderstatistics.jl",
               "mvnormal.jl",
               "mvnormalcanon.jl",
               "mvlognormal.jl",

src/univariate/orderstatistic.jl

@@ -0,0 +1,106 @@
+# Implementation based on chapters 2-4 of
+# Arnold, Barry C., Narayanaswamy Balakrishnan, and Haikady Navada Nagaraja.
+# A first course in order statistics. Society for Industrial and Applied Mathematics, 2008.
+
+"""
+    OrderStatistic{D<:UnivariateDistribution,S<:ValueSupport} <: UnivariateDistribution{S}
+
+The distribution of an order statistic from IID samples from a univariate distribution
+
+    OrderStatistic(dist::UnivariateDistribution, n::Int, i::Int; check_args::Bool=true)
+
+Construct the distribution of the `i`th order statistic from `n` independent samples from `dist`.
+
+The ``i``th order statistic of a sample is the ``i``th element of the sorted sample.
+For example, the 1st order statistic is the sample minimum, while the ``n``th order
+statistic is the sample maximum.
+
+If ``f`` is probability density (mass) function of `dist` with distribution function ``F``,
+then the probability density function ``g`` of the order statistic for continuous
+`dist` is
+```math
+g(x; n, i) = {n \\choose i} [F(x)]^{i-1} [1 - F(x)]^{n-i} f(x),
+```
+and the probability mass function ``g`` of the order statistic for discrete `dist` is
+```math
+g(x; n, i) = \\sum_{k=i}^n {n \\choose k} \\left( [F(x)]^k [1 - F(x)]^{n-k} - [F(x_-)]^k [1 - F(x_-)]^{n-k} \\right),
+```
+where ``x_-`` is the largest element in the support of `dist` less than ``x``.
+
+For the joint distribution of a subset of order statistics, use
+[`JointOrderStatistics`](@ref) instead.
+
+## Examples
+
+```julia
+OrderStatistic(Cauchy(), 10, 1)              # distribution of the sample minimum
+OrderStatistic(DiscreteUniform(10), 10, 10)  # distribution of the sample maximum
+OrderStatistic(Gamma(1, 1), 11, 5)           # distribution of the sample median
+```
+"""
+struct OrderStatistic{D<:UnivariateDistribution,S<:ValueSupport} <:
+       UnivariateDistribution{S}
+    dist::D
+    n::Int
+    i::Int
+    function OrderStatistic(
+        dist::UnivariateDistribution, n::Int, i::Int; check_args::Bool=true
+    )
+        @check_args(OrderStatistic, 1 ≤ i ≤ n)
+        return new{typeof(dist),value_support(typeof(dist))}(dist, n, i)
+    end
+end
+
+minimum(d::OrderStatistic) = minimum(d.dist)
+maximum(d::OrderStatistic) = maximum(d.dist)
+insupport(d::OrderStatistic, x::Real) = insupport(d.dist, x)
+
+params(d::OrderStatistic) = tuple(params(d.dist)..., d.n, d.i)
+partype(d::OrderStatistic) = partype(d.dist)
+Base.eltype(::Type{<:OrderStatistic{D}}) where {D} = Base.eltype(D)
+
+# distribution of the ith order statistic from an IID uniform distribution, with CDF Uᵢₙ(x)
+function _uniform_orderstatistic(d::OrderStatistic)
+    n = d.n
+    i = d.i
+    return Beta{Int}(i, n - i + 1)
+end
+
+function logpdf(d::OrderStatistic, x::Real)
+    b = _uniform_orderstatistic(d)
+    p = cdf(d.dist, x)
+    if value_support(typeof(d)) === Continuous
+        return logpdf(b, p) + logpdf(d.dist, x)
+    else
+        return logdiffcdf(b, p, p - pdf(d.dist, x))
+    end
+end
+
+for f in [:logcdf, :logccdf, :cdf, :ccdf]
+    @eval begin
+        function $f(d::OrderStatistic, x::Real)
+            b = _uniform_orderstatistic(d)
+            return $f(b, cdf(d.dist, x))
+        end
+    end
+end
+
+for f in [:quantile, :cquantile]
+    @eval begin
+        function $f(d::OrderStatistic, p::Real)
+            # since cdf is Fᵢₙ(x) = Uᵢₙ(Fₓ(x)), and Uᵢₙ is invertible and increasing, we
+            # have Fₓ(x) = Uᵢₙ⁻¹(Fᵢₙ(x)). then quantile function is
+            # Qᵢₙ(p) = inf{x: p ≤ Fᵢₙ(x)} = inf{x: Uᵢₙ⁻¹(p) ≤ Fₓ(x)} = Qₓ(Uᵢₙ⁻¹(p))
+            b = _uniform_orderstatistic(d)
+            return quantile(d.dist, $f(b, p))
+        end
+    end
+end
+
+function rand(rng::AbstractRNG, d::OrderStatistic)
+    # inverse transform sampling. Since quantile function is Qₓ(Uᵢₙ⁻¹(p)), we draw a random
+    # variable from Uᵢₙ and pass it through the quantile function of `d.dist`
+    T = eltype(d.dist)
+    b = _uniform_orderstatistic(d)
+    return T(quantile(d.dist, rand(rng, b)))
+end

src/univariates.jl

@@ -729,3 +729,5 @@ end
 for dname in continuous_distributions
     include(joinpath("univariate", "continuous", "$(dname).jl"))
 end
+
+include(joinpath("univariate", "orderstatistic.jl"))