[Fitting] Add confidence interval estimates through vcov (#311)

Co-authored-by: Oskar Laverny <[email protected]>

Author: Santymax98

Copulas.jl

@@ -9,7 +9,6 @@ Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 HCubature = "19dc6840-f33b-545b-b366-655c7e3ffd49"
-InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LambertW = "984bce1d-4616-540c-a9ee-88d1112d94c9"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"

Project.toml

@@ -39,7 +39,7 @@ plot(Ĉ)
 ### Full Model (with metadata)
 
 ```@example fitting_interface
-M = fit(CopulaModel, GumbelCopula, U; method=:default, summaries=true)
+M = fit(CopulaModel, GumbelCopula, U; method=:default)
 M
 ```
 
@@ -48,15 +48,15 @@ Returns a `CopulaModel` with:
 * `result` (the fitted copula), `n`, `ll` (log-likelihood),
 * `method`, `converged`, `iterations`, `elapsed_sec`,
 * `vcov` (if available),
-* `method_details` (a named tuple with method metadata and, if `summaries=true`, **pairwise summaries**: means, deviations, minima, and maxima of empirical τ/ρ/β/γ).
+* `method_details` (a named tuple with method metadata).
 
 ---
 
 ## Behavior & conventions (important)
 
 - fit operates on types, not on pre-constructed parameterized instances. Always pass a Copula or SklarDist *type* to `fit`, e.g. `fit(GumbelCopula, U)` or `fit(CopulaModel, SklarDist{ClaytonCopula,Tuple{Normal,LogNormal}}, X)`. If you already have a constructed instance `C0`, re-estimate its parameters by calling `fit(typeof(C0), U)`.
 
-- Default method selection: each family exposes the list of available fitting strategies via `_available_fitting_methods(CT)`. When `method = :default` the first element of that tuple is used. Example: `Copulas._available_fitting_methods(MyCopula)`.
+- Default method selection: each family exposes the list of available fitting strategies via `_available_fitting_methods(CT, d)`. When `method = :default` the first element of that tuple is used. Example: `Copulas._available_fitting_methods(MyCopula, d)`.
 
 - `CopulaModel` is the full result object returned by the fits performed via `Distributions.fit(::Type{CopulaModel}, ...)`. The light-weight shortcut `fit(MyCopula, U)` returns only a copula instance; use `fit(CopulaModel, ...)` to get diagnostics and metadata.
 
@@ -80,26 +80,7 @@ The `CopulaModel{CT} <: StatsBase.StatisticalModel` type stores the result and s
 | `hqc(M)`           | Hannan–Quinn criterion          |
 
 
-Quick access to the contained copula: `_copula_of(M)` (returns the copula even if `result` is a `SklarDist`).
-
-
-### Pairwise summaries and `method_details`
-
-When you request `summaries=true` (default) the returned `CopulaModel` contains extra pre-computed pairwise statistics inside `M.method_details`. Typical keys are:
-
-- `:tau_mean`, `:tau_sd`, `:tau_min`, `:tau_max`
-- `:rho_mean`, `:rho_sd`, `:rho_min`, `:rho_max`
-- `:beta_mean`, `:beta_sd`, `:beta_min`, `:beta_max`
-- `:gamma_mean`, `:gamma_sd`, `:gamma_min`, `:gamma_max`
-
-Access example:
-
-```@example fitting_interface
-M = fit(CopulaModel, GumbelCopula, U; summaries=true)
-M.method_details.tau_mean  # average pairwise Kendall's tau
-```
-
-If `summaries=false` these keys will be absent and `method_details` will be smaller.
+By default, the returned `CopulaModel` contains a lot of extra statistics, that you can see by printing the model in the REPL. 
 
 ### `vcov` and inference notes
 
@@ -143,7 +124,7 @@ plot(Ŝ.result)
 The names and availiability of fitting methods depends on the model. You can check what is available with the following internal call : 
 
 ```@example fitting_interface
-Copulas._available_fitting_methods(ClaytonCopula)
+Copulas._available_fitting_methods(ClaytonCopula, 3)
 ```
 
 The first method in the list is the one used by default. 
@@ -168,7 +149,7 @@ When you add a new copula family, implement the following so the generic `fit` f
 1. `_example(CT, d)` — return a representative instance (used to obtain default params and initial values).
 2. `_unbound_params(CT, d, params)` — transform the family `NamedTuple` parameters to an unconstrained `Vector{Float64}` used by optimizers.
 3. `_rebound_params(CT, d, α)` — invert `_unbound_params`, returning a `NamedTuple` suitable for `CT(d, ...)` construction.
-4. `_available_fitting_methods(::Type{<:YourCopula})` — declare supported methods (examples:  `:mle, :itau, :irho, :ibeta, ...`).
+4. `_available_fitting_methods(::Type{<:YourCopula}, d::Int)` — declare supported methods (examples:  `:mle, :itau, :irho, :ibeta, ...`).
 5. `_fit(::Type{<:YourCopula}, U, ::Val{:mle})` (and other `Val{}` methods) — implement the method and return `(fitted_copula, meta::NamedTuple)`; include keys such as `:θ̂`, `:optimizer`, `:converged`, `:iterations` and optionally `:vcov`.
 
 Place this checklist and a minimal `_fit` skeleton in `docs/src/manual/developer_fitting.md` where contributors can copy/paste and adapt.

docs/src/manual/fitting_interface.md

@@ -155,7 +155,7 @@ function _rebound_params(CT::Type{<:ArchimaxCopula{2, <:Generator, <:Tail}}, d,
     NamedTuple{all_names}(all_vals)
 end
 
-_available_fitting_methods(::Type{<:ArchimaxCopula}) = (:mle,)
+_available_fitting_methods(::Type{<:ArchimaxCopula}, d) = (:mle,)
 
 # 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]))
@@ -177,7 +177,7 @@ end
 
 # --- log-PDF stable ---
 function Distributions._logpdf(C::ArchimaxCopula{2, TG, TT}, u) where {TG, TT}
-    T = promote_type(Float64, eltype(u))
+    T = typeof(A(C.tail, one(ϕ(C.gen, one(eltype(u))))/2))
     @assert length(u) == 2
     u1, u2 = u
     (0.0 < u1 ≤ 1.0 && 0.0 < u2 ≤ 1.0) || return T(-Inf)

src/ArchimaxCopula.jl

@@ -191,11 +191,11 @@ _example(::Type{<:ArchimedeanCopula{d,<:FrailtyGenerator} where {d}}, d) = throw
 _unbound_params(CT::Type{<:ArchimedeanCopula}, d, θ) = _unbound_params(generatorof(CT), d, θ)
 _rebound_params(CT::Type{<:ArchimedeanCopula}, d, α) = _rebound_params(generatorof(CT), d, α)
 
-_available_fitting_methods(::Type{ArchimedeanCopula}) = (:gnz2011,)
-_available_fitting_methods(::Type{<:ArchimedeanCopula{d,GT} where {d,GT<:Generator}}) = (:mle,)
-_available_fitting_methods(::Type{<:ArchimedeanCopula{d,GT} where {d,GT<:UnivariateGenerator}}) = (:mle, :itau, :irho, :ibeta)
-_available_fitting_methods(::Type{<:ArchimedeanCopula{d,<:WilliamsonGenerator{d2, TX}} where {d,d2, TX}}) = Tuple{}() # No fitting method. 
-_available_fitting_methods(::Type{<:ArchimedeanCopula{d,<:WilliamsonGenerator{d2, <:Distributions.DiscreteNonParametric}} where {d,d2}}) = (:gnz2011,)
+_available_fitting_methods(::Type{ArchimedeanCopula}, d) = (:gnz2011,)
+_available_fitting_methods(::Type{<:ArchimedeanCopula{d,GT} where {d,GT<:Generator}}, d) = (:mle,)
+_available_fitting_methods(::Type{<:ArchimedeanCopula{d,GT} where {d,GT<:UnivariateGenerator}}, d) = (:mle, :itau, :irho, :ibeta)
+_available_fitting_methods(::Type{<:ArchimedeanCopula{d,<:WilliamsonGenerator{d2, TX}} where {d,d2, TX}}, d) = Tuple{}() # No fitting method. 
+_available_fitting_methods(::Type{<:ArchimedeanCopula{d,<:WilliamsonGenerator{d2, <:Distributions.DiscreteNonParametric}} where {d,d2}}, d) = (:gnz2011,)
 
 
 function _fit(::Union{Type{ArchimedeanCopula},Type{<:ArchimedeanCopula{d,<:WilliamsonGenerator{d2, <:Distributions.DiscreteNonParametric}} where {d,d2}}}, U, ::Val{:gnz2011})
@@ -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,27 +226,26 @@ 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)
     d = size(U,1)
     GT = generatorof(CT)
     lo, hi = _θ_bounds(GT, d)
-    θ₀ = [(lo+hi)/2]
+    θ₀ = [1.0]
     if start isa Real 
         θ₀[1] = start
     elseif start ∈ (:itau, :irho)
-        try 
-            θ₀[1] = only(Distributions.params(_fit(CT, U, Val{start}())[1]))
-        catch e
-        end
+        θ₀[1] = _fit(CT, U, Val{start}())[2].θ̂[1]
+    end
+    if θ₀[1] <= lo || θ₀[1] >= hi
+        θ₀[1] = Distributions.params(_example(CT, d))[1]
     end
-    θ₀[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/ArchimedeanCopula.jl

@@ -52,43 +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=10^4)
+    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
-    return H
+function ι(C::Copula{d}) where {d}
+    return Distributions.expectation(u -> -Distributions.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
@@ -128,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)

src/Copula.jl

@@ -2,7 +2,6 @@ module Copulas
 
     import Base
     import Random
-    import InteractiveUtils
     import SpecialFunctions
     import Roots
     import Distributions

src/Copulas.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}, d) = (:mle, :itau, :irho, :ibeta)

src/EllipticalCopulas/GaussianCopula.jl

@@ -97,4 +97,4 @@ function _rebound_params(::Type{<:TCopula}, d::Int, α::AbstractVector{T}) where
     return (; ν = ν, Σ = Σ)
 end
 
-_available_fitting_methods(::Type{<:TCopula}) = (:mle,)
+_available_fitting_methods(::Type{<:TCopula}, d) = (:mle,)

src/EllipticalCopulas/TCopula.jl

@@ -137,9 +137,9 @@ _example(CT::Type{<:ExtremeValueCopula}, d) = CT(d; _rebound_params(CT, d, fill(
 _unbound_params(CT::Type{<:ExtremeValueCopula}, d, θ) = _unbound_params(tailof(CT), d, θ)
 _rebound_params(CT::Type{<:ExtremeValueCopula}, d, α) = _rebound_params(tailof(CT), d, α)
 
-_available_fitting_methods(::Type{ExtremeValueCopula}) = (:ols, :cfg, :pickands)
-_available_fitting_methods(CT::Type{<:ExtremeValueCopula}) = (:mle,)
-_available_fitting_methods(CT::Type{<:ExtremeValueCopula{2,GT} where {GT<:UnivariateTail2}}) =  (:mle, :itau, :irho, :ibeta, :iupper)
+_available_fitting_methods(::Type{ExtremeValueCopula}, d) = (:ols, :cfg, :pickands)
+_available_fitting_methods(CT::Type{<:ExtremeValueCopula}, d) = (:mle,)
+_available_fitting_methods(CT::Type{<:ExtremeValueCopula{2,GT} where {GT<:UnivariateTail2}}, d) =  (:mle, :itau, :irho, :ibeta, :iupper)
 
 # Fitting empírico (OLS, CFG, Pickands):
 function _fit(::Type{ExtremeValueCopula}, U, method::Union{Val{:ols}, Val{:cfg}, Val{:pickands}}; 
@@ -152,25 +152,28 @@ function _fit(CT::Type{<:ExtremeValueCopula{d, GT} where {d, GT<:UnivariateTail2
         m isa Val{:irho} ? ρ⁻¹(CT,  StatsBase.corspearman(U')[1,2]) : 
                            β⁻¹(CT,  corblomqvist(U')[1,2])
     θ = clamp(θ, _θ_bounds(tailof(CT), 2)...)
-    return CT(2, θ), (; θ̂=θ)
+    return CT(2, θ), (; θ̂=(θ=θ,))
 end
 function _fit(CT::Type{<:ExtremeValueCopula{d, GT} where {d, GT<:UnivariateTail2}}, U, ::Val{:iupper})
     θ = clamp(λᵤ⁻¹(CT, λᵤ(U)), _θ_bounds(tailof(CT), 2)...)
-    return CT(2, θ), (; θ̂=θ)
+    return CT(2, θ), (; θ̂=(θ=θ,))
 end
 
 function _fit(CT::Type{<:ExtremeValueCopula{d, GT} where {d, GT<:UnivariateTail2}}, U, ::Val{:mle}; start::Union{Symbol,Real}=:itau, xtol::Real=1e-8)
     d = size(U,1)
     TT = tailof(CT)
     lo, hi = _θ_bounds(TT, d)
-    θ0 = start isa Real ? start : 
-         start ∈ (:itau, :irho, :ibeta, :iupper) ? _fit(CT, U, Val{start}())[2].θ̂ : 
-         only(Distributions.params(_example(CT, d)))
-    θ0 = clamp(θ0, lo, hi)
+    θ0_val = if start isa Real
+        start
+    else
+        initial_params = start ∈ (:itau, :irho, :ibeta, :iupper) ? _fit(CT, U, Val{start}())[2].θ̂ : only(Distributions.params(_example(CT, d)))
+        initial_params.θ
+    end
+    θ0_clamped = clamp(θ0_val, lo, hi)
     f(θ) = -Distributions.loglikelihood(CT(d, θ[1]), U)
-    res = Optim.optimize(f, lo, hi, [θ0], Optim.Fminbox(Optim.LBFGS()), autodiff = :forward)
+    res = Optim.optimize(f, lo, hi, [θ0_clamped], Optim.Fminbox(Optim.LBFGS()), autodiff = :forward)
     θ̂ = Optim.minimizer(res)[1]
-    return CT(d, θ̂), (; θ̂=θ̂, optimizer=:GradientDescent,
+    return CT(d, θ̂), (; θ̂=(;θ=θ̂), optimizer=:GradientDescent,
                         xtol=xtol, converged=Optim.converged(res), 
                         iterations=Optim.iterations(res))
 end