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

src/ArchimedeanCopula.jl

@@ -70,12 +70,25 @@ ArchimedeanCopula{D,TG}(d::Int, args...; kwargs...) where {D, TG} = ArchimedeanC
 
 Distributions.params(C::ArchimedeanCopula) = Distributions.params(C.G) # by default the parameter is the generator's parameters. 
 
-_cdf(C::ArchimedeanCopula, u) = ϕ(C.G, sum(ϕ⁻¹.(C.G, u)))
+function _cdf(C::ArchimedeanCopula, u) 
+    v = zero(eltype(u))
+    for uᵢ in u
+        v += ϕ⁻¹(C.G, uᵢ)
+    end
+    return ϕ(C.G, v)
+end
 function Distributions._logpdf(C::ArchimedeanCopula{d,TG}, u) where {d,TG}
-    if !all(0 .< u .< 1)
-        return eltype(u)(-Inf)
+    T = eltype(u)
+    for uᵢ in u
+        0 < uᵢ < 1 || return T(-Inf)
+    end
+    v = zero(eltype(u))
+    p = one(eltype(u))
+    for uᵢ in u
+        v += ϕ⁻¹(C.G, uᵢ)
+        p *= ϕ⁻¹⁽¹⁾(C.G, uᵢ)
     end
-    return log(max(ϕ⁽ᵏ⁾(C.G, d, sum(ϕ⁻¹.(C.G, u))) * prod(ϕ⁻¹⁽¹⁾.(C.G, u)), 0))
+    return log(max(ϕ⁽ᵏ⁾(C.G, d, v) * p, 0))
 end
 function Distributions._rand!(rng::Distributions.AbstractRNG, C::ArchimedeanCopula{d, TG}, x::AbstractVector{T}) where {T<:Real, d, TG}
     # By default, we use the Williamson sampling.
@@ -167,11 +180,19 @@ end
 function DistortionFromCop(C::ArchimedeanCopula, js::NTuple{p,Int}, uⱼₛ::NTuple{p,Float64}, i::Int) where {p}
     @assert length(js) == length(uⱼₛ)
     T = eltype(uⱼₛ)
-    sJ = sum(ϕ⁻¹.(C.G, uⱼₛ))
-    return ArchimedeanDistortion(C.G, p, float(sJ), float(T(ϕ⁽ᵏ⁾(C.G, p, sJ))))
+    sJ = zero(T)
+    for uⱼ in uⱼₛ
+        sJ += ϕ⁻¹(C.G, uⱼ)
+    end
+    rJ = ϕ⁽ᵏ⁾(C.G, p, sJ)
+    return ArchimedeanDistortion(C.G, p, float(sJ), float(rJ))
 end
 function ConditionalCopula(C::ArchimedeanCopula{D, TG}, ::NTuple{p,Int}, uⱼₛ::NTuple{p,Float64}) where {D, TG, p}
-    return ArchimedeanCopula(D - p, TiltedGenerator(C.G, p, sum(ϕ⁻¹.(C.G, uⱼₛ))))
+    sJ = zero(eltype(uⱼₛ))
+    for uⱼ in uⱼₛ
+        sJ += ϕ⁻¹(C.G, uⱼ)
+    end
+    return ArchimedeanCopula(D - p, TiltedGenerator(C.G, p, sJ))
 end
 SubsetCopula(C::ArchimedeanCopula{d,TG}, dims::NTuple{p, Int}) where {d,TG,p} = ArchimedeanCopula(length(dims), C.G)
 

src/Conditioning.jl

@@ -68,9 +68,11 @@ function Distributions.quantile(d::Distortion, α::Real)
     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: 
+# 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))
+Distributions.logpdf(d::Distortion, t::Real) = log(Distributions.pdf(d, t))
+Distributions.pdf(d::Distortion, t::Real) = exp(Distributions.logpdf(d, t))
 
 """
     DistortionFromCop{TC,p,T} <: Distortion
@@ -113,6 +115,21 @@ struct DistortionFromCop{TC,p}<: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
 
@@ -157,6 +174,45 @@ 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
 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).
+# Using copula calculus: f_{U_I|U_J}(u_I|u_J) = c(u)/den, and
+# f_{U_{i}|U_J}(u_i|u_J) = ∂^{p+1}C/∂(J,i) / den.
+function Distributions._logpdf(CC::ConditionalCopula{d,D,p,TDs}, v::AbstractVector{<:Real}) where {d,D,p,TDs}
+    T = promote_type(eltype(v), Float64)
+
+    # Support: 
+    CC.den <= 0 && return -Inf
+    for vₖ in v
+        0 < vₖ < 1 || return T(-Inf)
+    end
+
+    # 1) Map v → u_I via the stored distortions (non-sequential conditioning on J only)
+    uI = Distributions.quantile.(CC.distortions, v)
+    # 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: log c(u) - log den
+    logc_full = Distributions.logpdf(CC.C, u)
+    logden    = log(CC.den)
+    # 4) Sum of log conditional marginal numerators: ∑ log ∂^{p+1} C / ∂(J,i)
+    #    Evaluate each at vector with only (J fixed, i at u_i; others at 1)
+    sum_log_num = 0.0
+    for idx in 1:d
+        i = CC.is[idx]
+        ui = uI[idx]
+        z = _assemble(D, (i,), CC.js, (ui,), CC.uⱼₛ)
+        # Mixed partial of order p+1 with respect to (J..., i)
+        num = _der(u -> Distributions.cdf(CC.C, u), z, (CC.js..., i))
+        (num <= 0 || !isfinite(num)) && return T(-Inf)
+        sum_log_num += log(num)
+    end
+
+    # log c_{I|J}(v) = log c(u) + (d-1) log den − ∑ log num_i
+    return logc_full + (d - 1) * logden - sum_log_num
+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

src/Copula.jl

@@ -76,17 +76,39 @@ end
 function β(U::AbstractMatrix)
     # Assumes psuedo-data given. β multivariate (Hofert–Mächler–McNeil, ec. (7))
     d, n = size(U)
-    count = sum(j -> all(U[:, j] .<= 0.5) || all(U[:, j] .> 0.5), 1:n)
+    count = 0
+    @inbounds for j in 1:n
+        all_le = true; all_gt = true
+        for i in 1:d
+            v = U[i, j]
+            all_le &= (v <= 0.5)
+            all_gt &= (v > 0.5)
+            if !(all_le || all_gt)
+                break
+            end
+        end
+        count += (all_le || all_gt)
+    end
     h_d = 2.0^(d-1) / (2.0^(d-1) - 1.0)
     return h_d * (count/n - 2.0^(1-d))
 end
 function τ(U::AbstractMatrix)
     # Sample version of multivariate Kendall's tau for pseudo-data
     d, n = size(U)
     comp = 0
-    @inbounds for j in 2:n, i in 1:j-1
-        uᵢ = @view U[:, i]; uⱼ = @view U[:, j]
-        comp += (all(uᵢ .<= uⱼ) || all(uᵢ .>= uⱼ))
+    @inbounds for j in 2:n
+        for i in 1:j-1
+            le = true; ge = true
+            for k in 1:d
+                vij = U[k, i]; vjj = U[k, j]
+                le &= (vij <= vjj)
+                ge &= (vij >= vjj)
+                if !(le || ge)
+                    break
+                end
+            end
+            comp += (le || ge)
+        end
     end
     pc = comp / (n*(n-1)/2)
     return (2.0^d * pc - 2.0) / (2.0^d - 2.0)

src/EllipticalCopula.jl

@@ -50,22 +50,38 @@ abstract type EllipticalCopula{d,MT} <: Copula{d} end
 Base.eltype(C::CT) where CT<:EllipticalCopula = Base.eltype(N(CT)(C.Σ))
 function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, x::AbstractVector{T}) where {T<:Real, CT <: EllipticalCopula}
     Random.rand!(rng,N(CT)(C.Σ),x)
-    x .= clamp.(Distributions.cdf.(U(CT),x),0,1)
+    U₁ = U(CT)
+    @inbounds for i in eachindex(x)
+        x[i] = clamp(T(Distributions.cdf(U₁, x[i])), 0, 1)
+    end
     return x
 end
 function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, A::DenseMatrix{T}) where {T<:Real, CT<:EllipticalCopula}
     # More efficient version that precomputes stuff:
     n = N(CT)(C.Σ)
     u = U(CT)
     Random.rand!(rng,n,A)
-    A .= clamp.(Distributions.cdf.(u,A),0,1)
+    @inbounds for j in axes(A, 2), i in axes(A, 1)
+        A[i, j] = clamp(T(Distributions.cdf(u, A[i, j])), 0, 1)
+    end
     return A
 end
 function Distributions._logpdf(C::CT, u) where {CT <: EllipticalCopula}
     d = length(C)
     (u==zeros(d) || u==ones(d)) && return Inf 
-    x = StatsBase.quantile.(U(CT),u)
-    return Distributions.logpdf(N(CT)(C.Σ),x) - sum(Distributions.logpdf.(U(CT),x))
+    TΣ = eltype(C.Σ)
+    U₁ = U(CT)
+    # quantiles
+    x = Vector{TΣ}(undef, d)
+    @inbounds for i in 1:d
+        x[i] = StatsBase.quantile(U₁, u[i])
+    end
+    # sum of univariate logpdfs
+    s = zero(TΣ)
+    @inbounds for i in 1:d
+        s += Distributions.logpdf(U₁, x[i])
+    end
+    return Distributions.logpdf(N(CT)(C.Σ),x) - s
 end
 function make_cor!(Σ)
     # Verify that Σ is a correlation matrix, otherwise make it so : 

src/EllipticalCopulas/GaussianCopula.jl

@@ -52,7 +52,7 @@ struct GaussianCopula{d,MT} <: EllipticalCopula{d,MT}
             return IndependentCopula(size(Σ,1))
         end
         make_cor!(Σ)
-        N(GaussianCopula)(Σ)
+        Distributions.MvNormal(Σ) # only here for the checks... but is that really necessary ? 
         return new{size(Σ,1),typeof(Σ)}(Σ)
     end
 end
@@ -80,17 +80,16 @@ GaussianCopula{D, MT}(d::Int, ρ::Real) where {D, MT} = GaussianCopula(d, ρ)
 U(::Type{T}) where T<: GaussianCopula = Distributions.Normal()
 N(::Type{T}) where T<: GaussianCopula = Distributions.MvNormal
 function _cdf(C::CT,u) where {CT<:GaussianCopula}
-    x = StatsBase.quantile.(Distributions.Normal(), u)
     d = length(C)
-    return MvNormalCDF.mvnormcdf(C.Σ, fill(-Inf, d), x)[1]
+    return MvNormalCDF.mvnormcdf(C.Σ, fill(-Inf, d), StatsFuns.norminvcdf.(u))[1]
 end
 
 function rosenblatt(C::GaussianCopula, u::AbstractMatrix{<:Real})
-    return Distributions.cdf.(Distributions.Normal(), inv(LinearAlgebra.cholesky(C.Σ).L) * Distributions.quantile.(Distributions.Normal(), u))
+    return StatsFuns.normcdf.(inv(LinearAlgebra.cholesky(C.Σ).L) * StatsFuns.norminvcdf.(u))
 end
 
 function inverse_rosenblatt(C::GaussianCopula, s::AbstractMatrix{<:Real})
-    return Distributions.cdf.(Distributions.Normal(), LinearAlgebra.cholesky(C.Σ).L * Distributions.quantile.(Distributions.Normal(), s))
+    return StatsFuns.normcdf.(LinearAlgebra.cholesky(C.Σ).L * StatsFuns.norminvcdf.(s))
 end
 
 # Kendall tau of bivariate gaussian:
@@ -103,7 +102,7 @@ function DistortionFromCop(C::GaussianCopula{D,MT}, js::NTuple{p,Int}, uⱼₛ::
     ist = Tuple(setdiff(1:D, js))
     @assert i in ist
     J = collect(js)
-    zⱼ = Distributions.quantile.(Distributions.Normal(), collect(uⱼₛ))
+    zⱼ = StatsFuns.norminvcdf.(collect(uⱼₛ))
     if length(J) == 1 # if we condition on only one variable
         μz = C.Σ[i, J[1]] * zⱼ[1]
         σz = sqrt(1 - C.Σ[i, J[1]]^2)
@@ -137,7 +136,7 @@ function _rebound_params(::Type{<:GaussianCopula}, d::Int, α::AbstractVector{T}
     return (; Σ = _rebound_corr_params(d, α))
 end
 function _fit(CT::Type{<:GaussianCopula}, u, ::Val{:mle})
-    dd = Distributions.fit(N(CT), StatsBase.quantile.(U(CT),u))
+    dd = Distributions.fit(Distributions.MvNormal, StatsFuns.norminvcdf.(u))
     Σ = Matrix(dd.Σ)
     return GaussianCopula(Σ), (; θ̂ = (; Σ = Σ))
 end

src/ExtremeValueCopula.jl

@@ -47,7 +47,14 @@ ExtremeValueCopula{d,TT}(args...; kwargs...) where {d, TT} = ExtremeValueCopula(
 ExtremeValueCopula{D,TT}(d::Int, args...; kwargs...) where {D, TT} = ExtremeValueCopula{d,TT}(args...; kwargs...)
 (CT::Type{<:ExtremeValueCopula{2, <:Tail}})(d::Int, args...; kwargs...) = ExtremeValueCopula(2, tailof(CT)(args...; kwargs...))
 
-_cdf(C::ExtremeValueCopula{d, TT}, u) where {d, TT} = exp(-ℓ(C.tail, .- log.(u)))
+function _cdf(C::ExtremeValueCopula{d, TT}, u) where {d, TT}
+    d == length(u) || throw(ArgumentError("Dimension mismatch"))
+    z = similar(u)
+    @inbounds for i in 1:d
+        z[i] = -log(u[i])
+    end
+    return exp(-ℓ(C.tail, z))
+end
 Distributions.params(C::ExtremeValueCopula) = Distributions.params(C.tail)
 
 #### Restriction to bivariate cases of the following methods: 

src/Generator.jl

@@ -361,7 +361,7 @@ function ϕ⁽ᵏ⁾(G::WilliamsonGenerator{TX, d}, k::Int, t) where {d, TX<:Dis
     @inbounds for j in lastindex(r):-1:firstindex(r)
         rⱼ = r[j]; wⱼ = w[j]
         t ≥ rⱼ && break
-        zpow = (d == k+1) ? one(t) : (1 - t / rⱼ)^(d - 1 - k)
+        zpow = (d == k+1) ? one(Tt) : (1 - t / rⱼ)^(d - 1 - k)
         S += wⱼ * zpow / rⱼ^k
     end
     return S * (isodd(k) ? -1 : 1) * Base.factorial(d - 1) / Base.factorial(d - 1 - k)

src/Generator/ClaytonGenerator.jl

@@ -76,8 +76,11 @@ function Distributions._logpdf(C::ClaytonCopula{d,TG}, u) where {d,TG<:ClaytonGe
 
     θ = C.G.θ
     # Compute the sum of transformed variables
-    S1 = sum(t ^ (-θ) for t in u)
-    S2 = sum(log(t) for t in u)
+    S1 = S2 = zero(eltype(u))
+    for t in u
+        S1 += t^(-θ)
+        S2 += log(t)
+    end
     # Compute the log of the density according to the explicit formula for Clayton copula
     # See McNeil & Neslehova (2009), eq. (13)
     S1==d-1 && return eltype(u)(-Inf)

src/Generator/GumbelGenerator.jl

@@ -78,15 +78,16 @@ end
 
 function _cdf(C::ArchimedeanCopula{d,G}, u) where {d, G<:GumbelGenerator}
     θ = C.G.θ
-    lx = log.(.-log.(u))
-    return 1 - LogExpFunctions.cexpexp(LogExpFunctions.logsumexp(θ .* lx) ./ θ)
+    return 1 - LogExpFunctions.cexpexp(LogExpFunctions.logsumexp(θ .* log.( .- log.(u))) / θ)
 end
 function Distributions._logpdf(C::ArchimedeanCopula{2,GumbelGenerator{TF}}, u) where {TF}
     T = promote_type(TF, eltype(u))
-    !all(0 .< u .<= 1) && return T(-Inf) # if not in range return -Inf
+    u₁, u₂ = u 
+    (0 < u₁ <= 1) || return T(-Inf)
+    (0 < u₂ <= 1) || return T(-Inf)
 
     θ = C.G.θ
-    x₁, x₂ = -log(u[1]), -log(u[2])
+    x₁, x₂ = -log(u₁), -log(u₂)
     lx₁, lx₂ = log(x₁), log(x₂)
     A = LogExpFunctions.logaddexp(θ * lx₁, θ * lx₂)
     B = exp(A/θ)
@@ -105,29 +106,43 @@ end
 
 function Distributions._logpdf(C::ArchimedeanCopula{d,GumbelGenerator{TF}}, u) where {d,TF}
     T = promote_type(TF, eltype(u))
-    !all(0 .< u .<= 1) && return T(-Inf)
-
+    for uᵢ in u
+        (0 < uᵢ <= 1) || return T(-Inf)
+    end
+
     θ = C.G.θ
     α = 1 / θ
 
     # Step 1. Compute x_i = -log(u_i)
-    x = -log.(u)
-    lx = log.(x)
-
+    x = zeros(T, d)
+    θlx = zeros(T, d)
+    slx = zero(T)
+    sx = zero(T)
+    for (i, uᵢ) in enumerate(u)
+        xᵢ = -log(uᵢ)
+        lxᵢ = log(xᵢ)
+
+        x[i] = xᵢ
+        θlx[i] = lxᵢ * θ
+
+        sx += xᵢ
+        slx += lxᵢ
+    end
+
     # Step 2. Stable log-sum-exp for log(S) = log(sum(x_i^θ))
-    logS = LogExpFunctions.logsumexp(θ .* lx)
+    logt = LogExpFunctions.logsumexp(θlx)
 
     # Step 3. Compute log(φ⁽ᵈ⁾(S)) directly in log-domain
-    logt = logS
     tα = exp(α * logt)
     ntα = -tα
     logφ = -tα  # since log(φ(t)) = -exp(α * logt)
 
     # Compute the combinatorial sum in double precision
-    s = 0.0
+    s = zero(T)
     for j in 1:d
+
         term1 = α^j * Combinatorics.stirlings1(d, j, true)
-        inner_sum = 0.0
+        inner_sum = zero(T)
         for k in 1:j
             inner_sum += Combinatorics.stirlings2(j, k) * ntα^k
         end
@@ -141,7 +156,7 @@ function Distributions._logpdf(C::ArchimedeanCopula{d,GumbelGenerator{TF}}, u) w
     logφd = logφ - d * logt + log(abs(s))
 
     # Step 4. log|(φ⁻¹)'(u_i)| = log(θ) + (θ - 1)*log(-log(u_i)) - log(u_i)
-    sum_log_invderiv = d * log(θ) + (θ - 1)*sum(lx) - sum(log.(u))
+    sum_log_invderiv = d * log(θ) + (θ - 1) * slx + sx
 
     # Step 5. Combine
     return T(logφd + sum_log_invderiv)