first shot for add vcov

src/ArchimedeanCopula.jl

@@ -215,7 +215,7 @@ function _fit(CT::Type{<:ArchimedeanCopula{d, GT} where {d, GT<:UnivariateGenera
     θs = map(v -> invf(GT, clamp(v, -1, 1)), upper_triangle_flat)
 
     θ = clamp(Statistics.mean(θs), _θ_bounds(GT, d)...)
-    return CT(d, θ), (; θ̂=θ, eps)
+    return CT(d, θ), (; θ̂=(θ=θ,))
 end
 function _fit(CT::Type{<:ArchimedeanCopula{d, GT} where {d, GT<:UnivariateGenerator}}, U, ::Val{:ibeta})
     d    = size(U,1); δ = 1e-8; GT = generatorof(CT)
@@ -226,7 +226,7 @@ function _fit(CT::Type{<:ArchimedeanCopula{d, GT} where {d, GT<:UnivariateGenera
     βmin, βmax = fβ(a0), fβ(b0)
     if βmin > βmax; βmin, βmax = βmax, βmin; end
     θ = βobs ≤ βmin ? a0 : βobs ≥ βmax ? b0 : Roots.find_zero(θ -> fβ(θ)-βobs, (a0,b0), Roots.Brent(); xatol=1e-8, rtol=0)
-    return CT(d,θ), (; θ̂=θ)
+    return CT(d,θ), (; θ̂=(θ=θ,))
 end
 
 function _fit(CT::Type{<:ArchimedeanCopula{d, GT} where {d, GT<:UnivariateGenerator}}, U, ::Val{:mle}; start::Union{Symbol,Real}=:itau, xtol::Real=1e-8)
@@ -245,8 +245,8 @@ function _fit(CT::Type{<:ArchimedeanCopula{d, GT} where {d, GT<:UnivariateGenera
     θ₀[1] = clamp(θ₀[1], lo, hi)
     f(θ) = -Distributions.loglikelihood(CT(d, θ[1]), U)
     res = Optim.optimize(f, lo, hi,  θ₀, Optim.Fminbox(Optim.LBFGS()), autodiff = :forward)
-    θ̂     = Optim.minimizer(res)[1]
-    return CT(d, θ̂), (; θ̂=θ̂, optimizer=:GradientDescent,
+    θ     = Optim.minimizer(res)[1]
+    return CT(d, θ), (; θ̂=(θ=θ,), optimizer=Optim.summary(res),
                         xtol=xtol, converged=Optim.converged(res), 
                         iterations=Optim.iterations(res))
 end

src/Copula.jl

@@ -52,45 +52,15 @@ function β(C::Copula{d}) where {d}
     Cbar0 = Distributions.cdf(SurvivalCopula(C, Tuple(1:d)), u)
     return (2.0^(d-1) * C0 + Cbar0 - 1) / (2^(d-1) - 1)
 end
-function γ(C::Copula{d}; nmc::Int=100_000, rng::Random.AbstractRNG=Random.MersenneTwister(123)) where {d}
-    d ≥ 2 || throw(ArgumentError("γ(C) requires d≥2"))
-    if d == 2
-        f(t) = Distributions.cdf(C, [t, t]) + Distributions.cdf(C, [t, 1 - t])
-        I, _ = QuadGK.quadgk(f, 0.0, 1.0; rtol=sqrt(eps()))
-        return -2 + 4I
-    end
-    @inline _A(u)    = (minimum(u) + max(sum(u) - d + 1, 0.0)) / 2
-    @inline _Abar(u) = (1 - maximum(u) + max(1 - sum(u), 0.0)) / 2
-    @inline invfac(k::Integer) = exp(-SpecialFunctions.logfactorial(k))
-    s = 0.0
-    @inbounds for i in 0:d
-        s += (isodd(i) ? -1.0 : 1.0) * binomial(d, i) * invfac(i + 1)
-    end
-    a_d = 1/(d + 1) + 0.5*invfac(d + 1) + 0.5*s
-    b_d = 2/3 + 4.0^(1 - d) / 3
-    U = rand(rng, C, nmc)
-    m = 0.0
-    @inbounds for j in 1:nmc
-        u = @view U[:, j]
-        m += _A(u) + _Abar(u)
-    end
-    m /= nmc
-    return (m - a_d) / (b_d - a_d)
+function γ(C::Copula{d}) where {d}
+    _integrand(u) = (1 + minimum(u) - maximum(u) + max(abs(sum(u) - d/2) - (d - 2)/2, 0.0)) / 2
+    I = Distributions.expectation(_integrand, C; nsamples=1e5)
+    a = 1/(d+1) + 1/factorial(d+1)   # independence
+    b = (2 + 4.0^(1-d)) / 3          # comonotonicity
+    return (I - a) / (b - a)
 end
-
-function ι(C::Copula{d}; nmc::Int=100_000, rng::Random.AbstractRNG=Random.MersenneTwister(123)) where {d}
-    U = rand(rng, C, nmc)
-    s = 0.0
-    @inbounds for j in 1:nmc
-        u  = @view U[:, j]
-        lp = Distributions.logpdf(C, u)
-        isfinite(lp) || throw(DomainError(lp, "logpdf(C,u) non-finite."))
-        s -= lp
-    end
-    H = s / nmc
-    t = clamp(2H, -700.0, 0.0)
-    r = sqrt(max(0.0, 1 - exp(t)))
-    return (H = H, I = -H, r = r)
+function ι(C::Copula{d}) where {d}
+    return Distributions.expectation(u -> -logpdf(C, u), C; nsamples=10^4)
 end
 function λₗ(C::Copula{d}; ε::Float64 = 1e-10) where {d} 
     g(e) = Distributions.cdf(C, fill(e, d)) / e
@@ -130,51 +100,24 @@ function ρ(U::AbstractMatrix)
     return h * (2.0^d * μ - 1.0)
 end
 function γ(U::AbstractMatrix)
-    # Assumes pseudo-data given. Multivariate Gini’s gamma (Behboodian–Dolati–Úbeda, 2007)
     d, n = size(U)
-    if d == 2
-    # Schechtman–Yitzhaki symmetric Gini over ranks (copular invariant)
-        r1 = StatsBase.tiedrank(@view U[1, :])
-        r2 = StatsBase.tiedrank(@view U[2, :])
-        m  = n
-        h  = m + 1
-        acc = 0.0
-        @inbounds @simd for k in 1:m
-            acc += abs(r1[k] + r2[k] - h) - abs(r1[k] - r2[k])
-        end
-        return 2*acc / (m*h)
-    else
-        @inline _A(u)    = (minimum(u) + max(sum(u) - d + 1, 0.0)) / 2
-        @inline _Abar(u) = (1 - maximum(u) + max(1 - sum(u), 0.0)) / 2
-        invfac(k::Integer) = exp(-SpecialFunctions.logfactorial(k))
-        s = 0.0
-        binomf(d,i) = exp(SpecialFunctions.loggamma(d+1) - SpecialFunctions.loggamma(i+1) - SpecialFunctions.loggamma(d-i+1))
-        @inbounds for i in 0:d
-            s += (isodd(i) ? -1.0 : 1.0) * binomf(d,i) * invfac(i + 1)
-        end
-        a_d = 1/(d + 1) + 0.5*invfac(d + 1) + 0.5*s
-        b_d = 2/3 + 4.0^(1 - d) / 3
-        m = 0.0
-        @inbounds for j in 1:n
-            u = @view U[:, j]
-            m += _A(u) + _Abar(u)
-        end
-        m /= n
-        return (m - a_d) / (b_d - a_d)
+    I = zero(eltype(U))
+    for j in 1:n
+        u = U[:,j]
+        I += (1 + minimum(u) - maximum(u) + max(abs(sum(u) - d/2) - (d - 2)/2, 0.0)) / 2
     end
+    I /= n
+    a = 1/(d+1) + 1/factorial(d+1)
+    b = (2 + 4.0^(1-d)) / 3
+    return (I - a) / (b - a)
 end
 function _λ(U::AbstractMatrix; t::Symbol=:upper, p::Union{Nothing,Real}=nothing)
     # Assumes pseudo-data given. Multivariate tail’s lambda (Schmidt, R. & Stadtmüller, U. 2006)
-    d, m = size(U)
-    m ≥ 4 || throw(ArgumentError("At least 4 observations are required"))
-    p === nothing && (p = 1/sqrt(m))
+    p === nothing && (p = 1/sqrt(size(U, 2)))
     (0 < p < 1) || throw(ArgumentError("p must be in (0,1)"))
-    V = t === :upper ? (1 .- Float64.(U)) : Float64.(U)
-    cnt = 0
-    @inbounds @views for j in 1:m
-        cnt += all(V[:, j] .<= p)   # vista sin copiar gracias a @views
-    end
-    return clamp(cnt / (p*m), 0.0, 1.0)
+    in_tail = t=== :upper ? Base.Fix2(>=, 1-p) : Base.Fix2(<=, p)
+    prob = Statistics.mean(all(in_tail, U, dims=1))
+    return clamp(prob/p, 0.0, 1.0)
 end
 λₗ(U::AbstractMatrix; p::Union{Nothing,Real}=nothing) = _λ(U; t=:lower, p=p)
 λᵤ(U::AbstractMatrix; p::Union{Nothing,Real}=nothing) = _λ(U; t=:upper, p=p)
@@ -280,9 +223,7 @@ function ι(U::AbstractMatrix; k::Int=5, p::Real=Inf, leafsize::Int=32)
         logcd = d*log(2*SpecialFunctions.gamma(1 + 1/p)) - SpecialFunctions.loggamma(1 + d/p)
         H += logcd + (d / n) * sum(log.(ρ))
     end
-    t = clamp(2H, -700.0, 0.0)
-    r = sqrt(max(0.0, 1 - exp(t)))
-    return (H = H, I = -H, r = r)
+    return H
 end
 
 # Measure function. 

src/EllipticalCopulas/GaussianCopula.jl

@@ -82,7 +82,7 @@ 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, m=10_0000d)[1]
+    return MvNormalCDF.mvnormcdf(C.Σ, fill(-Inf, d), x)[1]
 end
 
 function rosenblatt(C::GaussianCopula, u::AbstractMatrix{<:Real})
@@ -139,6 +139,6 @@ end
 function _fit(CT::Type{<:GaussianCopula}, u, ::Val{:mle})
     dd = Distributions.fit(N(CT), StatsBase.quantile.(U(CT),u))
     Σ = Matrix(dd.Σ)
-    return GaussianCopula(Σ), (;)
+    return GaussianCopula(Σ), (; θ̂ = (; Σ = Σ))
 end
-_available_fitting_methods(::Type{<:GaussianCopula}) = (:itau, :irho, :ibeta, :mle)
+_available_fitting_methods(::Type{<:GaussianCopula}) = (:mle, :itau, :irho, :ibeta)