From 74a559c688dc51c4a029ef808afd127af919823f Mon Sep 17 00:00:00 2001 From: santymax98 Date: Fri, 28 Aug 2026 13:10:29 -0300 Subject: [PATCH] Add automatic copula family selection --- docs/src/manual/fitting_interface.md | 53 +++++++ src/Copulas.jl | 4 +- src/Fitting.jl | 221 +++++++++++++++++++++++++++ src/show.jl | 36 +++++ test/FittingSelectionTest.jl | 145 ++++++++++++++++++ test/runtests.jl | 1 + 6 files changed, 459 insertions(+), 1 deletion(-) create mode 100644 test/FittingSelectionTest.jl diff --git a/docs/src/manual/fitting_interface.md b/docs/src/manual/fitting_interface.md index 623469c59..65c1de5de 100644 --- a/docs/src/manual/fitting_interface.md +++ b/docs/src/manual/fitting_interface.md @@ -45,6 +45,59 @@ Returns a [`CopulaModel`](@ref) with: - `vcov` (if available), - `method_details` (a `NamedTuple` with method-specific metadata). +## Automatic copula-family selection + +When the copula family is unknown, `CopulaModel` can select it automatically +from a collection of candidate families: + +```@example fitting_interface +Ctrue = ClaytonCopula(2, 4.0) +Usel = rand(Ctrue, 1_000) + +Msel = fit( + CopulaModel, + Copula, + Usel; + candidates=(ClaytonCopula, GumbelCopula, FrankCopula), + criterion=:bic, + vcov=false, +) +Msel +``` + +The available information criteria are: + +- `:bic` — Bayesian information criterion, +- `:aic` — Akaike information criterion, +- `:aicc` — finite-sample corrected AIC, +- `:hqc` — Hannan–Quinn criterion. + +With `criterion=:default`, BIC is used. + +Candidate fits are compared using the requested criterion, and the family with +the smallest eligible finite value is selected. The winning family is then +fitted once more using the inference options requested by the user. + +The complete comparison can be inspected with [`selectiontable`](@ref): + +```@example fitting_interface +selectiontable(Msel) +``` + +Each row stores the candidate family, fitting status and method, +log-likelihood, number of parameters, and all four information criteria. +Candidates that fail to fit can be skipped with `on_error=:skip` (the default) +or propagated immediately with `on_error=:throw`. + +Two built-in candidate collections are available: + +- `candidates=:default` uses a conservative set of commonly used parametric families; +- `candidates=:all` considers the broader built-in parametric repertoire compatible with the data dimension. + +An explicit tuple of families is recommended when the scientific problem +already restricts the plausible candidate set. + + --- ## Behavior & conventions (important) diff --git a/src/Copulas.jl b/src/Copulas.jl index d5d160d21..f6e24512c 100644 --- a/src/Copulas.jl +++ b/src/Copulas.jl @@ -148,6 +148,7 @@ module Copulas 𝒲, TiltedGenerator, EmpiricalGenerator, + Copula, SklarDist, # SklarDist to make multivariate models AMHCopula, # And a bunch of copulas. ArchimedeanCopula, @@ -198,6 +199,7 @@ module Copulas BernsteinCopula, BetaCopula, CheckerboardCopula, - CopulaModel + CopulaModel, + selectiontable end diff --git a/src/Fitting.jl b/src/Fitting.jl index d97b2b3dc..a3f833972 100644 --- a/src/Fitting.jl +++ b/src/Fitting.jl @@ -653,3 +653,224 @@ function StatsBase.predict(M::CopulaModel; newdata=nothing, what=:cdf, nsim=0) what === :pdf ? (newdata === nothing ? throw(ArgumentError("`newdata` required for `:pdf`")) : Distributions.pdf(C, newdata)) : throw(ArgumentError("`what` must be one of :simulate, :cdf, or :pdf. Got `$what`.")) end + +############################################################################### +##### Automatic copula-family selection +############################################################################### + +""" + _available_selection_criteria() + +Return the information criteria available for automatic copula selection. +The first entry is used when `criterion=:default`. +""" +_available_selection_criteria() = (:bic, :aic, :aicc, :hqc) + +struct CopulaSelectionTable{T,R<:AbstractVector{T}} <: AbstractVector{T} + rows::R + criterion::Symbol + selected_family +end + +Base.IndexStyle(::Type{<:CopulaSelectionTable}) = IndexLinear() +Base.size(table::CopulaSelectionTable) = size(table.rows) +Base.getindex(table::CopulaSelectionTable, i::Int) = table.rows[i] +Base.sort(table::CopulaSelectionTable; kwargs...) = + CopulaSelectionTable(sort(table.rows; kwargs...), table.criterion, table.selected_family) + +""" + _default_copula_candidates(d) + +Return the built-in parametric copula families considered by automatic +selection in dimension `d` when `candidates=:default`. +""" +function _default_copula_candidates(d::Integer) + common = ( + GaussianCopula, + TCopula, + AMHCopula, + ClaytonCopula, + FrankCopula, + GumbelCopula, + JoeCopula, + GalambosCopula, + HuslerReissCopula, + LogCopula, + ) + d == 2 && return (common..., PlackettCopula) + return common +end + +""" + _all_copula_candidates(d) + +Return the broad built-in parametric copula repertoire considered by automatic +selection when `candidates=:all`. +""" +function _all_copula_candidates(d::Integer) + common = ( + IndependentCopula, + GaussianCopula, + TCopula, + AMHCopula, + ClaytonCopula, + FrankCopula, + GumbelCopula, + GumbelBarnettCopula, + InvGaussianCopula, + JoeCopula, + BB1Copula, + BB2Copula, + BB3Copula, + BB6Copula, + BB7Copula, + BB8Copula, + BB9Copula, + BB10Copula, + RafteryCopula, + GalambosCopula, + HuslerReissCopula, + LogCopula, + CuadrasAugeCopula, + MixedCopula, + MOCopula, + TawnCopula, + ) + d == 2 && return ( + common..., + PlackettCopula, + FGMCopula, + BB4Copula, + BB5Copula, + AsymGalambosCopula, + AsymLogCopula, + AsymMixedCopula, + BC2Copula, + ) + return common +end + +function _selection_candidates(candidates, d::Integer) + candidates === :default && return _default_copula_candidates(d) + candidates === :all && return _all_copula_candidates(d) + return (candidates isa Type || candidates isa UnionAll) ? (candidates,) : Tuple(candidates) +end + +""" + selectiontable(model::CopulaModel) + +Return the candidate-comparison table stored in an automatically selected +copula model. +""" +function selectiontable(M::CopulaModel) + get(M.method_details, :selection, false) || + throw(ArgumentError("The model was not produced by automatic copula selection.")) + return CopulaSelectionTable( + M.method_details.selection_table, + M.method_details.criterion, + M.method_details.selected_family, + ) +end + +""" + fit(CopulaModel, Copula, U; candidates=:default, criterion=:default, method=:default, kwargs...) + +Fit candidate copula families to the `d x n` pseudo-observation matrix `U`, +select the family minimizing the requested information criterion, and return +the selected model. + +Candidate fits used only for comparison are performed with `vcov=false`. The +winning family is fitted once with the inference options requested by the user. +""" +function Distributions.fit(::Type{CopulaModel}, ::Type{Copula}, U; + candidates=:default, criterion::Symbol=:default, method::Symbol=:default, + on_error::Symbol=:skip, require_convergence::Bool=true, + quick_fit::Bool=false, derived_measures::Bool=true, vcov::Bool=true, + vcov_method=nothing, kwargs...) + d, _ = size(U) + available_criteria = _available_selection_criteria() + criterion = criterion === :default ? first(available_criteria) : criterion + criterion in available_criteria || + throw(ArgumentError("Criterion '$criterion' is not available. Available: $(join(available_criteria, ", ")).")) + on_error in (:skip, :throw) || throw(ArgumentError("`on_error` must be either `:skip` or `:throw`.")) + + candidate_types = _selection_candidates(candidates, d) + + rows = NamedTuple[] + best_type = nothing + best_method = nothing + best_score = Inf + best_index = 0 + selection_start = time() + + for CT in candidate_types + CT <: Copula || throw(ArgumentError("Candidate `$CT` is not a subtype of `Copula`.")) + CT === Copula && throw(ArgumentError("`Copula` cannot itself appear inside `candidates`.")) + try + M = Distributions.fit(CopulaModel, CT, U; method=method, + quick_fit=false, derived_measures=false, vcov=false, kwargs...) + aic_value = StatsBase.aic(M) + aicc_value = aicc(M) + bic_value = StatsBase.bic(M) + hqc_value = hqc(M) + score = + criterion === :bic ? bic_value : + criterion === :aic ? aic_value : + criterion === :aicc ? aicc_value : hqc_value + status = + !isfinite(M.ll) || !isfinite(score) ? :nonfinite : + require_convergence && !M.converged ? :not_converged : + :ok + if status === :nonfinite + aic_value = isfinite(aic_value) ? aic_value : Inf + aicc_value = isfinite(aicc_value) ? aicc_value : Inf + bic_value = isfinite(bic_value) ? bic_value : Inf + hqc_value = isfinite(hqc_value) ? hqc_value : Inf + score = Inf + end + push!(rows, (candidate=CT, status=status, method=M.method, + converged=M.converged, nparams=StatsBase.dof(M), + loglikelihood=M.ll, aic=aic_value, aicc=aicc_value, + bic=bic_value, hqc=hqc_value, + criterion_value=score, elapsed_sec=M.elapsed_sec, + error=nothing)) + if status === :ok && score < best_score + best_type = CT + best_method = M.method + best_score = score + best_index = length(rows) + end + catch err + on_error === :throw && rethrow() + push!(rows, (candidate=CT, status=:failed, method=method, + converged=false, nparams=0, loglikelihood=NaN, + aic=Inf, aicc=Inf, bic=Inf, hqc=Inf, + criterion_value=Inf, elapsed_sec=NaN, + error=sprint(showerror, err))) + end + end + + best_type === nothing && throw(ErrorException("No candidate copula produced an eligible finite fit.")) + selected = Distributions.fit(CopulaModel, best_type, U; method=best_method, + quick_fit=quick_fit, derived_measures=derived_measures, vcov=vcov, + vcov_method=vcov_method, kwargs...) + quick_fit && return selected + + elapsed_sec = time() - selection_start + return CopulaModel(selected.result, selected.n, selected.ll, selected.method; + vcov=selected.vcov, + converged=selected.converged, + iterations=selected.iterations, + elapsed_sec=elapsed_sec, + method_details=(; selected.method_details..., + selection=true, + criterion=criterion, + candidates=candidate_types, + requested_method=method, + selection_options=(; on_error, require_convergence, kwargs...), + selection_table=rows, + selected_family=best_type, + selected_index=best_index, + selected_score=best_score, + selection_elapsed_sec=elapsed_sec)) +end diff --git a/src/show.jl b/src/show.jl index 1fcf38b7a..17f44d6f2 100644 --- a/src/show.jl +++ b/src/show.jl @@ -64,6 +64,31 @@ Pretty p-value formatting: show very small values as inequalities. """ _pstr(p) = p < 1e-16 ? "<1e-16" : Printf.@sprintf("%.4g", p) +function Base.show(io::IO, ::MIME"text/plain", table::CopulaSelectionTable) + criterion_label = uppercase(String(table.criterion)) + println(io, "Copula model selection (criterion: ", criterion_label, ")") + _hr(io) + Printf.@printf(io, "%-2s %-24s %-14s %-10s %12s %12s\n", + "", "Family", "Status", "Method", "LogLik", criterion_label) + + for row in table + selected = row.candidate === table.selected_family ? "-" : "" + family = replace(string(row.candidate), "Copulas." => "") + loglik = isfinite(row.loglikelihood) ? + Printf.@sprintf("%.3f", row.loglikelihood) : + string(row.loglikelihood) + value = getproperty(row, table.criterion) + criterion_value = isfinite(value) ? Printf.@sprintf("%.3f", value) : "-" + Printf.@printf(io, "%-2s %-24s %-14s %-10s %12s %12s\n", + selected, family, String(row.status), String(row.method), + loglik, criterion_value) + end + + _hr(io) + println(io, "- selected model") +end + + """ Key-value aligned printing for header lines. """ @@ -183,6 +208,17 @@ function Base.show(io::IO, M::CopulaModel) end _kv(io, "Number of observations", Printf.@sprintf("%d", StatsBase.nobs(M))) + md = M.method_details + if get(md, :selection, false) + _section(io, "Model selection") + _kv(io, "Criterion", uppercase(String(md.criterion))) + _kv(io, "Candidate families", string(length(md.candidates))) + _kv(io, "Selected family", string(md.selected_family)) + excluded = count(row -> row.status !== :ok, md.selection_table) + excluded > 0 && _kv(io, "Excluded or failed", string(excluded)) + end + + _section(io, "Fit metrics") ll = M.ll ll0 = get(M.method_details, :null_ll, NaN) diff --git a/test/FittingSelectionTest.jl b/test/FittingSelectionTest.jl new file mode 100644 index 000000000..85eec7542 --- /dev/null +++ b/test/FittingSelectionTest.jl @@ -0,0 +1,145 @@ +@testset "Automatic copula-family selection" begin + U = rand(rng, ClaytonCopula(2, 6.0), 300) + + candidates = ( + IndependentCopula, + ClaytonCopula, + ) + + @testset "BIC selection" begin + M = fit( + CopulaModel, + Copula, + U; + candidates=candidates, + criterion=:bic, + vcov=false, + ) + + @test M isa CopulaModel + @test M.method_details.selection === true + @test M.method_details.criterion === :bic + @test M.method_details.candidates == candidates + @test M.method_details.selected_family === ClaytonCopula + @test M.method_details.selected_index in eachindex(candidates) + @test isfinite(M.method_details.selected_score) + + table = selectiontable(M) + + @test table isa AbstractVector + @test length(table) == length(candidates) + @test table.criterion === :bic + @test table.selected_family === ClaytonCopula + + @test [row.candidate for row in table] == collect(candidates) + @test all(row -> row.status === :ok, table) + @test all(row -> isfinite(row.loglikelihood), table) + @test all(row -> isfinite(row.bic), table) + + selected_row = only(filter( + row -> row.candidate === ClaytonCopula, + table, + )) + + @test selected_row.bic == minimum(row.bic for row in table) + @test M.method_details.selected_score == selected_row.bic + end + + @testset "Selection table display" begin + M = fit( + CopulaModel, + Copula, + U; + candidates=candidates, + vcov=false, + ) + + table = selectiontable(M) + + io = IOBuffer() + show(io, MIME("text/plain"), table) + printed = String(take!(io)) + + @test occursin("Copula model selection", printed) + @test occursin("BIC", printed) + @test occursin("IndependentCopula", printed) + @test occursin("ClaytonCopula", printed) + @test occursin("selected model", printed) + + io = IOBuffer() + show(io, MIME("text/plain"), M) + printed = String(take!(io)) + + @test occursin("Model selection", printed) + @test occursin("Criterion:", printed) + @test occursin("Selected family:", printed) + end + + @testset "Information criteria" begin + for criterion in (:aic, :aicc, :hqc) + M = fit( + CopulaModel, + Copula, + U; + candidates=candidates, + criterion=criterion, + vcov=false, + ) + + table = selectiontable(M) + + @test M.method_details.criterion === criterion + @test table.criterion === criterion + + eligible = filter(row -> row.status === :ok, table) + values = getproperty.(eligible, criterion) + + @test M.method_details.selected_score == minimum(values) + end + end + + @testset "API errors" begin + fitted = fit( + CopulaModel, + ClaytonCopula, + U; + vcov=false, + ) + + @test_throws ArgumentError selectiontable(fitted) + + @test_throws ArgumentError fit( + CopulaModel, + Copula, + U; + candidates=candidates, + criterion=:invalid, + vcov=false, + ) + + @test_throws ArgumentError fit( + CopulaModel, + Copula, + U; + candidates=candidates, + on_error=:invalid, + vcov=false, + ) + + @test_throws ArgumentError fit( + CopulaModel, + Copula, + U; + candidates=(Normal,), + vcov=false, + ) + + @test_throws ArgumentError fit( + CopulaModel, + Copula, + U; + candidates=(Copula,), + vcov=false, + ) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 4679ba562..8c99eb5e7 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -16,6 +16,7 @@ testfiles = [ "EllipticalCopulas", "ExpectationMaximizationExt", "FittingTest", + "FittingSelectionTest", "MiscelaneousCopulas", "NatafTest", "SklarDist",