[New Copulas] d > 2 EV and Archimax

src/ArchimaxCopula.jl

@@ -40,10 +40,45 @@ ArchimaxCopula(d, gen::Generator, ::NoTail) = ArchimedeanCopula(d, gen)
 ArchimaxCopula(d, ::IndependentGenerator, tail::Tail) = ExtremeValueCopula(d, tail) 
 Distributions.params(C::ArchimaxCopula) = (_as_tuple(Distributions.params(C.gen))..., _as_tuple(Distributions.params(C.tail))...)
 
-# Fast conditional distortion binding (bivariate)
-DistortionFromCop(C::ArchimaxCopula{2}, js::NTuple{1,Int}, uⱼₛ::NTuple{1,Float64}, ::Int) = BivArchimaxDistortion(C.gen, C.tail, Int8(js[1]), float(uⱼₛ[1]))
+function _cdf(C::ArchimaxCopula{d}, u) where {d}
+    # Basic support checks
+    T = eltype(u)
+    any(iszero, u) && return T(0)
+    all(isone, u) && return T(1)
+    # Compute x_i = ϕ⁻¹(u_i), S = ∑ x_i, ω = x/S, then ϕ(S·A(ω))
+    x = ϕ⁻¹.(C.gen, u)
+    S = sum(x)
+    S == 0 && return T(1)
+    ω = ntuple(i -> x[i] / S, d)
+    return ϕ(C.gen, S * A(C.tail, ω))
+end
+function Distributions._logpdf(C::ArchimaxCopula{d, TG, TT}, u) where {d, TG, TT}
+    @inbounds for ui in u
+        (0.0 < ui < 1.0) || return -Inf
+    end
+    val = _der(v -> Distributions.cdf(C, v), collect(u), ntuple(identity, d))
+    return log(max(val, 0))
+end
+function Distributions._rand!(rng::Distributions.AbstractRNG, C::ArchimaxCopula{d, TG, TT}, X::AbstractMatrix{T}) where {T<:Real, d, TG, TT}
+    d == 2 && return @invoke Distributions._rand!(rng::Distributions.AbstractRNG, C::ArchimaxCopula{2, TG, TT}, X)
+    @assert size(X, 1) == d
+    U = rand(rng, Distributions.Uniform(), d, size(X, 2))
+    X .= inverse_rosenblatt(C, U)
+    return X
+end
+function Distributions._rand!(rng::Distributions.AbstractRNG, C::ArchimaxCopula{d, TG, TT}, x::AbstractVector{T}) where {T<:Real, d, TG, TT}
+    d == 2 && return @invoke Distributions._rand!(rng::Distributions.AbstractRNG, C::ArchimaxCopula{2, TG, TT}, x)
+    u = rand(rng, Distributions.Uniform(), d)
+    x .= inverse_rosenblatt(C, u)
+    return x
+end
 
-# --- CDF ---
+
+
+
+
+
+###### Special methods for the bivariate cases
 function _cdf(C::ArchimaxCopula{2}, u)
     u1, u2 = u
     (0.0 ≤ u1 ≤ 1.0 && 0.0 ≤ u2 ≤ 1.0) || return 0.0
@@ -57,8 +92,6 @@ function _cdf(C::ArchimaxCopula{2}, u)
     t = _safett(y / S)                 # protect t≈0,1
     return ϕ(C.gen, S * A(C.tail, t))
 end
-
-# --- log-PDF stable ---
 function Distributions._logpdf(C::ArchimaxCopula{2, TG, TT}, u) where {TG, TT}
     T = promote_type(Float64, eltype(u))
     @assert length(u) == 2
@@ -90,17 +123,6 @@ function Distributions._logpdf(C::ArchimaxCopula{2, TG, TT}, u) where {TG, TT}
     base > 0 || return T(-Inf)
     return T(log(φpp) + log(base))
 end
-
-# --- Kendall τ: τ = τ_A + (1 - τ_A) τ_ψ ---
-τ(C::ArchimaxCopula) = begin
-    τA = τ(ExtremeValueCopula(2, C.tail))
-    τψ = τ(C.gen)
-    τA + (1 - τA) * τψ
-end
-
-
-# Use the matrix sampler for better efficiency
-# (if not working, maybe uncomment the vetor version ?)
 function Distributions._rand!(rng::Distributions.AbstractRNG, C::ArchimaxCopula{2, TG, TT}, A::DenseMatrix{T}) where {T<:Real, TG, TT}
     evcop, frail = ExtremeValueCopula(2, C.tail), frailty(C.gen)
     Distributions._rand!(rng, evcop, A)
@@ -116,6 +138,12 @@ function Distributions._rand!(rng::Distributions.AbstractRNG, C::ArchimaxCopula{
     x[2] = ϕ(C.gen, -log(v2)/M)
     return x
 end
+DistortionFromCop(C::ArchimaxCopula{2}, js::NTuple{1,Int}, uⱼₛ::NTuple{1,Float64}, ::Int) = BivArchimaxDistortion(C.gen, C.tail, Int8(js[1]), float(uⱼₛ[1]))
+τ(C::ArchimaxCopula{2, TG, TT})  where {TG, TT} = begin
+    τA = τ(ExtremeValueCopula(2, C.tail))
+    τψ = τ(C.gen)
+    τA + (1 - τA) * τψ
+end
 
 
 

src/ExtremeValueCopula.jl

@@ -45,7 +45,29 @@ end
 _cdf(C::ExtremeValueCopula{d, TT}, u) where {d, TT} = exp(-ℓ(C.tail, .- log.(u)))
 Distributions.params(C::ExtremeValueCopula) = Distributions.params(C.tail)
 
-#### Restriction to bivariate cases of the following methods: 
+function Distributions._rand!(rng::Distributions.AbstractRNG, C::ExtremeValueCopula{d, TT}, X::AbstractMatrix{T}) where {T<:Real, d, TT}
+    @assert size(X, 1) == d
+    U = rand(rng,  d, size(X, 2))
+    X .= inverse_rosenblatt(C, U)
+    return X
+end
+function Distributions._rand!(rng::Distributions.AbstractRNG, C::ExtremeValueCopula{d, TT}, x::AbstractVector{T}) where {T<:Real, d, TT}
+    u = rand(rng, d)
+    x .= inverse_rosenblatt(C, u)
+    return x
+end
+function Distributions._logpdf(C::ExtremeValueCopula{d, TT}, u) where {d, TT}
+    # domain checks
+    @inbounds for ui in u
+        (0.0 < ui < 1.0) || return -Inf
+    end
+    # Compute mixed partial ∂^d/∂u₁…∂u_d of the cdf via nested ForwardDiff
+    val = _der(v -> Distributions.cdf(C, v), collect(u), 1:d)
+    return log(max(val, 0))
+end
+
+
+###### Restriction to bivariate cases of the following methods: 
 function Distributions._logpdf(C::ExtremeValueCopula{2, TT}, u) where {TT}
     u1, u2 = u
     (0.0 < u1 ≤ 1.0 && 0.0 < u2 ≤ 1.0) || return -Inf 
@@ -86,4 +108,4 @@ function Distributions._rand!(rng::Distributions.AbstractRNG, C::ExtremeValueCop
     x[2] = exp(log(w)*(1-z)/a)
     return x
 end
-DistortionFromCop(C::ExtremeValueCopula{2, TT}, js::NTuple{1,Int}, uⱼₛ::NTuple{1,Float64}, ::Int) where TT = BivEVDistortion(C.tail, Int8(js[1]), float(uⱼₛ[1]))
+DistortionFromCop(C::ExtremeValueCopula{d, TT}, js::NTuple{1,Int}, uⱼₛ::NTuple{1,Float64}, ::Int) where {d, TT} = BivEVDistortion(C.tail, Int8(js[1]), float(uⱼₛ[1]))

src/Tail.jl

@@ -33,21 +33,18 @@ References:
 abstract type Tail end
 Base.broadcastable(tail::Tail) = Ref(tail)
 
-####### Functions you need to overload: 
+####### Main interface for any dimensions: 
 _is_valid_in_dim(tail::Tail, d::Int) = throw(ArgumentError("Validity of the tail type $(typeof(tail)) must be supplied by overwriting the function _is_valid_in_dim(tail::Tail, d::Int)"))
 A(::Tail, ω::NTuple{d,<:Real}) where {d} = throw(ArgumentError("Implement A(Tail{$d}, ω) en el simplex Δ_{d-1}"))
-
-####### Rest of the interface you can overload if more efficient:
-needs_binary_search(::Tail) = false
-# \ell function
 function ℓ(tail::Tail, x)
     s = sum(x)
     return s == 0 ? zero(eltype(x)) : s * A(tail, ntuple(i->x[i]/s, length(x)))
 end
 
 
-# A more friendly interface for models that are only bivariate: 
+####### A more friendly interface for models that are only bivariate: 
 abstract type Tail2 <: Tail end
+needs_binary_search(::Tail2) = false
 _is_valid_in_dim(::Tail2, d::Int) = (d==2)
 A(tail::Tail2, t::NTuple{2, <:Real}) = A(tail, t[1])
 dA(tail::Tail2, t::Real) = ForwardDiff.derivative(z -> A(tail, z), t)

src/Tail/GalambosTail.jl

@@ -38,16 +38,20 @@ end
 const GalambosCopula{T} = ExtremeValueCopula{2, GalambosTail{T}}
 GalambosCopula(θ) =ExtremeValueCopula(2, GalambosTail(θ))
 Distributions.params(tail::GalambosTail) = (tail.θ,)
+_is_valid_in_dim(::GalambosTail, d::Int) = (d >= 2)
 
+function A(tail::GalambosTail, ω::NTuple{d,<:Real}) where {d}
+    θ = tail.θ
+    θ == 0 && return 1.0
+    isinf(θ) && return maximum(ω)
+    return -LogExpFunctions.expm1(-LogExpFunctions.logsumexp(-θ .* log.(ω))/θ)  # 1 - (∑ ω_i^{-θ})^{-1/θ}
+end
+
+#### Special bindings for dimension d == 2
 needs_binary_search(tail::GalambosTail) = (tail.θ > 19.5)
 function A(tail::GalambosTail, t::Real)
     tt = _safett(t)
-    θ  = tail.θ
-    if θ == 0
-        return 1.0
-    elseif isinf(θ)
-        return max(tt, 1-tt)
-    else
-        return -LogExpFunctions.expm1(-LogExpFunctions.logaddexp(-θ*log(tt), -θ*log(1-tt)) / θ)
-    end
-end
+    tail.θ == 0 && return 1.0
+    isinf(tail.θ) && return max(tt, 1-tt)
+    return -LogExpFunctions.expm1(-LogExpFunctions.logaddexp(-tail.θ*log(tt), -tail.θ*log(1-tt)) / tail.θ)
+end