[Feat] Add Distortions's pdf/logpdf, unroll a few loops.

Project.toml

@@ -1,6 +1,6 @@
 name = "Copulas"
 uuid = "ae264745-0b69-425e-9d9d-cf662c5eec93"
-version = "0.1.36"
+version = "0.1.37"
 authors = ["Oskar Laverny"]
 
 [deps]

src/Conditioning.jl

@@ -63,22 +63,25 @@ function Distributions.quantile(d::Distortion, α::Real)
     T = typeof(float(α))
     ϵ = eps(T)
     α < ϵ && return zero(T)
-    α > 1 - 2ϵ && return one(T)
-    lα = log(α)
-    f(u) = Distributions.logcdf(d, u) - lα
-    return Roots.find_zero(f, (ϵ, 1 - 2ϵ), Roots.Bisection(); xtol = sqrt(eps(T)))
-end
-# You have to implement one of these two: 
+    α > 1 - 2ϵ && return one(T)
+    lα = log(α)
+    f(u) = Distributions.logcdf(d, u) - lα
+    lo, hi = ϵ, one(T) - 2ϵ
+    flo, fhi = f(lo), f(hi)
+    flo >= zero(flo) && return lo
+    fhi <= zero(fhi) && return hi
+    return Roots.find_zero(f, (lo, hi), Roots.Bisection(); xtol = sqrt(eps(T)))
+end
+# You have to implement a cdf, and you can implement a pdf, either in log scaleor not: 
 Distributions.logcdf(d::Distortion, t::Real) = log(Distributions.cdf(d, t))
 Distributions.cdf(d::Distortion, t::Real) = exp(Distributions.logcdf(d, t))
-
-# These slow versions are given, but you should probably overrid them: 
 function Distributions.logpdf(d::Distortion, u::Real)
     (0.0 <= u <= 1.0) || return -Inf
     v = ForwardDiff.derivative(t -> Distributions.cdf(d, t), float(u))
     v <= 0 && return -Inf
     return log(v)
 end
+Distributions.pdf(d::Distortion, t::Real) = exp(Distributions.logpdf(d, t))
 
 """
     DistortionFromCop{TC,p,T} <: Distortion
@@ -119,6 +122,21 @@ struct DistortionFromCop{TC,p,T}<:Distortion
 end
 Distributions.cdf(d::DistortionFromCop, u::Real) = _partial_cdf(d.C, (d.i,), d.js, (u,), d.uⱼₛ) / d.den
 
+# Density on the uniform scale: f_{i|J}(u | u_J) = (∂^{p+1} C / ∂(J..., i))(u, u_J) / (∂^{p} C / ∂J)(1, u_J)
+function Distributions.logpdf(d::DistortionFromCop, u::Real)
+    # Support checks
+    (0 < u < 1) || return -Inf
+    d.den <= 0 && return -Inf
+
+    # Assemble the evaluation point with u at coordinate i, u_J fixed, others at 1
+    D = length(d.C)
+    z = _assemble(D, (d.i,), d.js, (float(u),), d.uⱼₛ)
+    # Mixed partial derivative of order p+1 w.r.t. (J..., i)
+    num = _der(u -> Distributions.cdf(d.C, u), z, (d.js..., d.i))
+    (num <= 0 || !isfinite(num)) && return -Inf
+    return log(num) - log(d.den)
+end
+
 """
     DistortedDist{Disto,Distrib} <: Distributions.UnivariateDistribution
 
@@ -132,6 +150,7 @@ struct DistortedDist{Disto, Distrib}<:Distributions.ContinuousUnivariateDistribu
     end
 end
 Distributions.cdf(D::DistortedDist, t::Real) = Distributions.cdf(D.D, Distributions.cdf(D.X, t))
+Distributions.logcdf(D::DistortedDist, t::Real) = Distributions.logcdf(D.D, Distributions.cdf(D.X, t))
 Distributions.quantile(D::DistortedDist, α::Real) = Distributions.quantile(D.X, Distributions.quantile(D.D, α))
 function Distributions.logpdf(D::DistortedDist, t::Real)
     u = Distributions.cdf(D.X, t)
@@ -149,22 +168,56 @@ struct ConditionalCopula{d, D, p, T, TDs}<:Copula{d}
     is::NTuple{d, Int}
     uⱼₛ::NTuple{p, T}
     den::T
+    logden::T
     distortions::TDs
     function ConditionalCopula(C::Copula{D}, js, uⱼₛ) where {D}
         jst, uⱼₛt = _process_tuples(Val{D}(), js, uⱼₛ)
-        ist = Tuple(setdiff(1:D, jst))
+        ist = Tuple(i for i in 1:D if i ∉ jst)
         p = length(jst)
         d = D - p
         distos = Tuple(DistortionFromCop(C, jst, uⱼₛt, i) for i in ist)
-        den = p==1 ? Distributions.pdf(subsetdims(C, jst), uⱼₛt[1]) :
-                     Distributions.pdf(subsetdims(C, jst), collect(uⱼₛt))
+        den = all(disto -> disto isa DistortionFromCop, distos) ? distos[1].den :
+              (p==1 ? Distributions.pdf(subsetdims(C, jst), uⱼₛt[1]) :
+                      Distributions.pdf(subsetdims(C, jst), collect(uⱼₛt)))
         T = promote_type(eltype(uⱼₛt), typeof(den))
-        return new{d, D, p, T, typeof(distos)}(C, jst, ist, NTuple{p,T}(uⱼₛt), T(den), distos)
+        denT = T(den)
+        return new{d, D, p, T, typeof(distos)}(
+            C, jst, ist, NTuple{p,T}(uⱼₛt), denT,
+            denT > zero(T) ? log(denT) : T(-Inf), distos
+        )
     end
 end
 function _cdf(CC::ConditionalCopula{d,D,p,T}, v::AbstractVector{<:Real}) where {d,D,p,T}
-    return _partial_cdf(CC.C, CC.is, CC.js, Distributions.quantile.(CC.distortions, v), CC.uⱼₛ) / CC.den
+    uI = ntuple(k -> Distributions.quantile(CC.distortions[k], v[k]), d)
+    return _partial_cdf(CC.C, CC.is, CC.js, uI, CC.uⱼₛ) / CC.den
+end
+
+# Density of the conditional copula on [0,1]^d.
+# For v ∈ (0,1)^d, let u_I = quantile(D_k, v_k) with D_k = H_{i_k|J}(·|u_J).
+# Then c_{I|J}(v) = f_{U_I|U_J}(u_I|u_J) / ∏_k f_{U_{i_k}|U_J}(u_{i_k}|u_J).
+# The Jacobian of v_k = D_k(u_k) contributes 1 / pdf(D_k, u_k).
+function Distributions._logpdf(CC::ConditionalCopula{d,D,p,T,TDs}, v::AbstractVector{<:Real}) where {d,D,p,T,TDs}
+    TR = promote_type(eltype(v), T)
+
+    # Support:
+    CC.den <= 0 && return TR(-Inf)
+    for vₖ in v
+        0 < vₖ < 1 || return TR(-Inf)
+    end
+
+    # 1) Map v → u_I via the stored distortions (non-sequential conditioning on J only)
+    uI = ntuple(k -> Distributions.quantile(CC.distortions[k], v[k]), d)
+    # 2) Full u vector at which to evaluate the base copula density
+    u = _assemble(D, CC.is, CC.js, uI, CC.uⱼₛ)
+    # 3) Joint conditional density on the original uniform scale
+    logdensity = TR(Distributions.logpdf(CC.C, u) - CC.logden)
+    # 4) Change variables from u_I to the conditional marginal scales v
+    for idx in 1:d
+        logdensity -= Distributions.logpdf(CC.distortions[idx], uI[idx])
+    end
+    return logdensity
 end
+
 # Sampling: sequential inverse-CDF using conditional distortions
 function Distributions._rand!(rng::Distributions.AbstractRNG, CC::ConditionalCopula{d, D, p}, x::AbstractVector{T}) where {T<:Real, d, D, p}
     # We want a sample from the COPULA of the conditional model. Let U be a
@@ -235,18 +288,30 @@ condition(C::Copula{D}, j, xⱼ) where D = condition(C, _process_tuples(Val{D}()
 # (StackOverflow). The downstream `DistortionFromCop`/`ConditionalCopula` still
 # store `Float64`, so non-`Float64` values are converted there — the conditioning
 # result is computed in `Float64` regardless of input precision.
+function _conditional_components(C::Copula, js, uⱼₛ, is)
+    CC = ConditionalCopula(C, js, uⱼₛ)
+    distortions = CC isa ConditionalCopula ? CC.distortions :
+                  Tuple(DistortionFromCop(C, js, uⱼₛ, i) for i in is)
+    return CC, distortions
+end
+
 function condition(C::Copula{D}, js::NTuple{p, Int}, uⱼₛ::NTuple{p, <:Real}) where {D, p}
-    margins = Tuple(DistortionFromCop(C, js, uⱼₛ, i) for i in setdiff(1:D, js))
-    p==D-1 && return margins[1]
-    return SklarDist(ConditionalCopula(C, js, uⱼₛ), margins)
+    is = Tuple(setdiff(1:D, js))
+    p==D-1 && return DistortionFromCop(C, js, uⱼₛ, is[1])
+    CC, distortions = _conditional_components(C, js, uⱼₛ, is)
+    return SklarDist(CC, distortions)
 end
 
 condition(C::SklarDist{<:Copula{D}}, j, xⱼ) where D = condition(C, _process_tuples(Val{D}(), j, xⱼ)...)
 function condition(X::SklarDist{<:Copula{D}, Tpl}, js::NTuple{p, Int}, xⱼₛ::NTuple{p, <:Real}) where {D, Tpl, p}
     uⱼₛ = Tuple(Distributions.cdf(X.m[j], xⱼ) for (j,xⱼ) in zip(js, xⱼₛ))
-    margins = Tuple(DistortionFromCop(X.C, js, uⱼₛ, i)(X.m[i]) for i in setdiff(1:D, js))
-    p==D-1 && return margins[1]
-    return SklarDist(ConditionalCopula(X.C, js, uⱼₛ), margins)
+    is = Tuple(setdiff(1:D, js))
+    if p == D - 1
+        return DistortionFromCop(X.C, js, uⱼₛ, is[1])(X.m[is[1]])
+    end
+    CC, distortions = _conditional_components(X.C, js, uⱼₛ, is)
+    margins = Tuple(distortions[k](X.m[is[k]]) for k in eachindex(is))
+    return SklarDist(CC, margins)
 end
 
 ###########################################################################
@@ -348,4 +413,4 @@ function inverse_rosenblatt(D::SklarDist, u::AbstractMatrix{<:Real})
         v[i,:] .= Distributions.quantile.(Mᵢ, v[i,:])
     end
     return v
-end
+end

src/EllipticalCopulas/GaussianCopula.jl

@@ -86,7 +86,9 @@ function _cdf(C::CT,u) where {CT<:GaussianCopula}
 end
 
 function rosenblatt(C::GaussianCopula, u::AbstractMatrix{<:Real})
-    return Distributions.cdf.(Distributions.Normal(), inv(LinearAlgebra.cholesky(C.Σ).L) * Distributions.quantile.(Distributions.Normal(), u))
+    L = LinearAlgebra.cholesky(C.Σ).L
+    z = Distributions.quantile.(Distributions.Normal(), u)
+    return Distributions.cdf.(Distributions.Normal(), L \ z)
 end
 
 function inverse_rosenblatt(C::GaussianCopula, s::AbstractMatrix{<:Real})
@@ -122,6 +124,23 @@ function ConditionalCopula(C::GaussianCopula{D,MT}, js::NTuple{p,Int}, uⱼₛ::
     return GaussianCopula(Σcond)
 end
 
+function _conditional_components(C::GaussianCopula{D,MT}, js::NTuple{p,Int},
+                                 uⱼₛ::NTuple{p,Float64}, is) where {D,MT,p}
+    J = collect(Int, js)
+    I = collect(Int, is)
+    Σ = C.Σ
+    F = LinearAlgebra.cholesky(LinearAlgebra.Symmetric(Σ[J, J]))
+    zJ = Distributions.quantile.(Distributions.Normal(), collect(uⱼₛ))
+    ΣIJ = Σ[I, J]
+    μ = ΣIJ * (F \ zJ)
+    Σcond = Σ[I, I] - ΣIJ * (F \ Σ[J, I])
+    distortions = ntuple(k -> begin
+        σ² = max(Σcond[k, k], zero(eltype(Σcond)))
+        GaussianDistortion(float(μ[k]), float(sqrt(σ²)))
+    end, length(is))
+    return GaussianCopula(Σcond), distortions
+end
+
 # Subsetting colocated
 SubsetCopula(C::GaussianCopula, dims::NTuple{p, Int}) where p = GaussianCopula(C.Σ[collect(dims),collect(dims)])
 
@@ -141,4 +160,4 @@ function _fit(CT::Type{<:GaussianCopula}, u, ::Val{:mle})
     Σ = Matrix(dd.Σ)
     return GaussianCopula(Σ), (; θ̂ = (; Σ = Σ))
 end
-_available_fitting_methods(::Type{<:GaussianCopula}, d) = (:mle, :itau, :irho, :ibeta)
+_available_fitting_methods(::Type{<:GaussianCopula}, d) = (:mle, :itau, :irho, :ibeta)

src/EllipticalCopulas/TCopula.jl

@@ -42,11 +42,69 @@ TCopula{D,df,MT}(d::Int, ν::Real, Σ::AbstractMatrix)  where {D,df,MT} = TCopul
 
 
 U(::Type{TCopula{d,df,MT}}) where {d,df,MT} = Distributions.TDist(df)
-N(::Type{TCopula{d,df,MT}}) where {d,df,MT} = function(Σ)
-    Distributions.MvTDist(df,Σ)
-end
-
-# Kendall tau of bivariate student: 
+N(::Type{TCopula{d,df,MT}}) where {d,df,MT} = function(Σ)
+    Distributions.MvTDist(df,Σ)
+end
+
+function _student_rosenblatt_cache(C::TCopula{d}) where d
+    Σ = C.Σ
+    return ntuple(d) do k
+        k == 1 && return nothing
+        J = 1:(k - 1)
+        F = LinearAlgebra.cholesky(LinearAlgebra.Symmetric(Σ[J, J]))
+        β = F \ Σ[J, k]
+        σ0² = max(Σ[k, k] - LinearAlgebra.dot(Σ[k, J], β), zero(eltype(Σ)))
+        return (; F, β, σ0 = sqrt(σ0²))
+    end
+end
+
+function rosenblatt(C::TCopula{d,ν}, u::AbstractMatrix{<:Real}) where {d,ν}
+    size(u, 1) == d || throw(ArgumentError("Dimension mismatch between copula and input matrix"))
+    Tu = Distributions.TDist(ν)
+    z = Distributions.quantile.(Tu, u)
+    v = similar(z)
+    v[1, :] .= u[1, :]
+    cache = _student_rosenblatt_cache(C)
+    @inbounds for k in 2:d
+        entry = cache[k]
+        Tcond = Distributions.TDist(ν + k - 1)
+        for col in axes(u, 2)
+            zJ = view(z, 1:(k - 1), col)
+            solved_zJ = entry.F \ zJ
+            μ = LinearAlgebra.dot(entry.β, zJ)
+            δ = LinearAlgebra.dot(zJ, solved_zJ)
+            σ = entry.σ0 * sqrt((ν + δ) / (ν + k - 1))
+            v[k, col] = Distributions.cdf(Tcond, (z[k, col] - μ) / σ)
+        end
+    end
+    return v
+end
+
+function inverse_rosenblatt(C::TCopula{d,ν}, s::AbstractMatrix{<:Real}) where {d,ν}
+    size(s, 1) == d || throw(ArgumentError("Dimension mismatch between copula and input matrix"))
+    Tu = Distributions.TDist(ν)
+    z = similar(s, float(promote_type(eltype(s), eltype(C.Σ))))
+    v = similar(z)
+    z[1, :] .= Distributions.quantile.(Tu, s[1, :])
+    v[1, :] .= s[1, :]
+    cache = _student_rosenblatt_cache(C)
+    @inbounds for k in 2:d
+        entry = cache[k]
+        Tcond = Distributions.TDist(ν + k - 1)
+        for col in axes(s, 2)
+            zJ = view(z, 1:(k - 1), col)
+            solved_zJ = entry.F \ zJ
+            μ = LinearAlgebra.dot(entry.β, zJ)
+            δ = LinearAlgebra.dot(zJ, solved_zJ)
+            σ = entry.σ0 * sqrt((ν + δ) / (ν + k - 1))
+            z[k, col] = μ + σ * Distributions.quantile(Tcond, s[k, col])
+            v[k, col] = Distributions.cdf(Tu, z[k, col])
+        end
+    end
+    return v
+end
+
+# Kendall tau of bivariate student:
 # Lindskog, F., McNeil, A., & Schmock, U. (2003). Kendall’s tau for elliptical distributions. In Credit risk: Measurement, evaluation and management (pp. 149-156). Heidelberg: Physica-Verlag HD.
 τ(C::TCopula{2,MT}) where MT = 2*asin(C.Σ[1,2])/π 
 
@@ -58,10 +116,11 @@ function DistortionFromCop(C::TCopula{D,ν,MT}, js::NTuple{p,Int}, uⱼₛ::NTup
     if length(Jv) == 1
         r = RiJ[1]; μz = r * zJ[1]; σ0² = 1 - r^2; δ = zJ[1]^2
     else
-        L = LinearAlgebra.cholesky(LinearAlgebra.Symmetric(ΣJJ))
-        μz = LinearAlgebra.dot(RiJ, (L' \ (L \ zJ)))
-        σ0² = 1 - LinearAlgebra.dot(RiJ, (L' \ (L \ RJi)))
-        y = L \ zJ; δ = LinearAlgebra.dot(y, y)
+        F = LinearAlgebra.cholesky(LinearAlgebra.Symmetric(ΣJJ))
+        solved_zJ = F \ zJ
+        μz = LinearAlgebra.dot(RiJ, solved_zJ)
+        σ0² = 1 - LinearAlgebra.dot(RiJ, F \ RJi)
+        δ = LinearAlgebra.dot(zJ, solved_zJ)
     end
     νp = ν + length(Jv); σz = sqrt(max(σ0², zero(σ0²))) * sqrt((ν + δ) / νp)
     return StudentDistortion(float(μz), float(σz), Int(ν), Int(νp))
@@ -78,6 +137,29 @@ function ConditionalCopula(C::TCopula{D,df,MT}, js, uⱼₛ) where {D,df,MT}
     R_cond = Matrix(Σcond ./ (σ * σ'))
     return TCopula(df + p, R_cond)
 end
+
+function _conditional_components(C::TCopula{D,ν,MT}, js::NTuple{p,Int},
+                                 uⱼₛ::NTuple{p,Float64}, is) where {D,ν,MT,p}
+    J = collect(Int, js)
+    I = collect(Int, is)
+    Σ = C.Σ
+    F = LinearAlgebra.cholesky(LinearAlgebra.Symmetric(Σ[J, J]))
+    zJ = Distributions.quantile.(Distributions.TDist(ν), collect(uⱼₛ))
+    solved_zJ = F \ zJ
+    ΣIJ = Σ[I, J]
+    μ = ΣIJ * solved_zJ
+    Σcond = Σ[I, I] - ΣIJ * (F \ Σ[J, I])
+    δ = LinearAlgebra.dot(zJ, solved_zJ)
+    νp = ν + p
+    scale = sqrt((ν + δ) / νp)
+    distortions = ntuple(k -> begin
+        σ² = max(Σcond[k, k], zero(eltype(Σcond)))
+        StudentDistortion(float(μ[k]), float(sqrt(σ²) * scale), Int(ν), Int(νp))
+    end, length(is))
+    σ = sqrt.(LinearAlgebra.diag(Σcond))
+    Rcond = Matrix(Σcond ./ (σ * σ'))
+    return TCopula(νp, Rcond), distortions
+end
 # Subsetting colocated
 SubsetCopula(C::TCopula{d,df,MT}, dims::NTuple{p, Int}) where {d,df,MT,p} = TCopula(df, C.Σ[collect(dims),collect(dims)])
 

src/Generator.jl

@@ -45,10 +45,49 @@ max_monotony(G::Generator) = throw("This generator does not have a defined max m
 ϕ⁽¹⁾(G::Generator, t) = ForwardDiff.derivative(x -> ϕ(G,x), t)
 ϕ⁻¹⁽¹⁾(G::Generator, t) = ForwardDiff.derivative(x -> ϕ⁻¹(G, x), t)
 ϕ⁽ᵏ⁾(G::Generator, k::Int, t) = taylor(ϕ(G), t, k)[end] * factorial(k)
-ϕ⁽ᵏ⁾⁻¹(G::Generator, k::Int, t; start_at=t) = try 
-    Roots.find_zero(x -> ϕ⁽ᵏ⁾(G, k, x) - t, start_at)
-catch
-    Roots.find_zero(x -> ϕ⁽ᵏ⁾(G, k, x) - t, (0,Inf))
+function ϕ⁽ᵏ⁾⁻¹(G::Generator, k::Int, t; start_at=t)
+    f(x) = ϕ⁽ᵏ⁾(G, k, x) - t
+    T = typeof(float(t))
+    lo, hi = eps(T), one(T)
+    flo, fhi = f(lo), f(hi)
+    iszero(flo) && return lo
+    iszero(fhi) && return hi
+
+    for _ in 1:64
+        signbit(flo) != signbit(fhi) &&
+            return Roots.find_zero(f, (lo, hi), Roots.Bisection())
+        hi *= 2
+        fhi = f(hi)
+        iszero(fhi) && return hi
+    end
+    throw(ArgumentError("Could not bracket the inverse generator derivative"))
+end
+
+# Stable evaluations of W(exp(logx)) and W₋₁(-exp(logx)). They avoid forming
+# arguments that overflow or underflow close to independence in generators
+# whose first-derivative inverses have a Lambert-W closed form.
+function _lambertw_exp(logx::T) where {T<:AbstractFloat}
+    isinf(logx) && return logx > 0 ? T(Inf) : zero(T)
+    logx <= log(floatmax(T)) && return LambertW.lambertw(exp(logx))
+
+    w = logx - log(logx)
+    for _ in 1:4
+        w -= (w + log(w) - logx) / (one(T) + inv(w))
+    end
+    return w
+end
+
+function _lambertwm1_negexp(logabsx::T) where {T<:AbstractFloat}
+    logabsx == T(-Inf) && return T(-Inf)
+    logabsx = min(logabsx, -one(T))
+    logabsx >= log(floatmin(T)) && return LambertW.lambertw(-exp(logabsx), -1)
+
+    target = -logabsx
+    y = target + log(target)
+    for _ in 1:4
+        y -= (y - log(y) - target) / (one(T) - inv(y))
+    end
+    return -y
 end