Copulas.jl 1.0: from a Distributions.jl experiment to a stable dependence-modeling toolkit

Copulas.jl has reached version 1.0. The documentation now describes the stable public API, the package’s current model bestiary, and the boundary between supported behavior and implementation details.

Copulas.jl 1.0 is now released.

The most important aspect of this 1.0 release is not a single new algorithm. The package has gained many models and operations, but the main milestone is more fundamental:

Copulas.jl now has a clear public contract that can be maintained under semantic versioning.

Copulas.jl started from a simple design idea: copulas are distribution functions, so in Julia they should behave like distributions. More than four years later, that idea still anchors the package, while the surrounding implementation and API have expanded substantially.

This post is therefore both a 1.0 announcement and a summary of the package’s evolution.

The original idea

When Copulas.jl was first announced in 2022, the goal was to provide a native Julia implementation of standard copula workflows while fitting naturally into the Distributions.jl ecosystem.

Two abstractions were already central:

• Copula, the common distribution type for dependence models;
• SklarDist, which combines a copula and arbitrary univariate margins into a multivariate distribution through Sklar’s theorem.

The attraction of that design was not only syntactic consistency. Once a copula or a Sklar distribution behaves like a normal Julia distribution, a lot of ecosystem interoperability comes almost for free.

Sampling is rand. Densities are pdf and logpdf. Distribution functions are cdf. Fitting builds on Distributions.fit. Likelihood-based downstream packages can consume the same objects instead of requiring a special copula adapter at every boundary.

That became visible very early. In 2022, a generic density implementation for Archimedean copulas made it possible to put a SklarDist directly inside a Turing model and use the normal likelihood interface. That was a small feature in terms of API surface, but it was a useful demonstration of the direction of the package: build the mathematical object once, and let Julia’s interfaces do the composition.

The same philosophy still drives 1.0.

The Archimedean detour that became a core strength

A lot of the early development focused on Archimedean copulas.

Named families such as Clayton, Frank, Gumbel, Joe, Ali–Mikhail–Haq, and inverse Gaussian were progressively added, but the more interesting question was always how much of their implementation could be shared.

That led to the generator-based architecture and, eventually, to Williamson-transform machinery.

By the 0.1.17 announcement in 2023, Copulas.jl could use the Williamson (d)-transform to obtain radial representations for generic Archimedean generators. This made it possible to provide generic sampling paths even when there was no family-specific stochastic representation available.

That machinery has kept evolving since then:

• more named and BB-family generators;
• frailty-based fast paths when complete monotonicity gives a simpler representation;
• numerical and exact inverse-Williamson constructions;
• real Williamson orders;
• exact beta-product reductions when changing dimension/order;
• removal of factorial overflow in high-dimensional derivative and radial calculations;
• and, in 1.0, an explicitly documented public extension protocol for user-defined generators.

The last point matters a lot for the meaning of the 1.0 release.

A downstream user can now define a subtype of Generator and implement the documented mathematical contract — in particular ϕ and max_monotony, together with the documented parameter interface when the generator is parameterized — and rely on that as supported public API. The many inverse/derivative methods, radial caches, fitting hooks, dispatch traits, reconstruction machinery, and optimization tricks used internally by Copulas.jl remain internal.

That is exactly the kind of distinction that was previously implicit and is now explicit.

From a bestiary to operations on dependence models

For a long time, adding copula families was the most visible kind of progress. But the bigger transformation of Copulas.jl happened when the package started adding operations that work across families.

The 0.1.31 era brought a major batch of them:

• subsetting and marginalization;
• conditioning;
• Rosenblatt transforms;
• inverse Rosenblatt transforms;
• richer conditional-distribution machinery;
• plotting recipes;
• BB Archimedean families;
• and the first generic Archimax API.

This changed what one could do with a copula object.

Instead of only asking “can I sample from this family?” or “can I evaluate its density?”, one could start composing workflows:

using Copulas, Distributions

C = ClaytonCopula(3, 2.0)

C23 = subsetdims(C, (2, 3))
Cgiven = condition(C, 2, 0.4)

u = rand(C, 100)
s = rosenblatt(C, u)
u2 = inverse_rosenblatt(C, s)

The exact implementation may differ dramatically between Gaussian, Student, Archimedean, extreme-value, empirical, or composed models, but the user-level vocabulary stays the same.

Conditionals in particular are just distributions. The stable interface is condition followed by ordinary operations such as cdf, pdf, quantile, or rand; family-specific h-functions and inverse h-functions can remain optimized implementation details.

This brings the package closer to a general toolkit for dependence modeling rather than a collection of family-specific implementations.

Models beyond the classical named families

The model bestiary also kept expanding.

By the time of 1.0, the package includes a broad collection of:

• Archimedean models, including Clayton, Frank, Gumbel, Joe, AMH, inverse Gaussian, and BB1–BB10;
• multivariate Gaussian and Student copulas;
• extreme-value families such as Logistic, Galambos, Hüsler–Reiss, extremal-t, Tawn, asymmetric and spectral constructions;
• empirical, empirical-beta, Bernstein, checkerboard, and empirical extreme-value models;
• bivariate Archimax constructions;
• nested Archimedean copulas;
• Liouville copulas;
• as well as independence, Fréchet bounds, survival copulas, subsets, and related constructions.

Some of these were previously announced; many were not.

The Discourse thread is actually a nice record of the project changing in public. In late 2023, someone asked when nested copulas would arrive. At the time, the honest answer was that they were planned but still required quite a bit of machinery.

They are now here.

NestedArchimedeanCopula supports hierarchical Archimedean constructions of arbitrary tree depth, with generator choices at the subtree level. The implementation reuses the package’s common generator and Taylor-derivative infrastructure instead of building a parallel mathematical stack just for nested models.

Version 1.0 closes that loop.

Extreme values grew up too

The last major announcement on the original Discourse thread before this one was 0.1.24, which introduced bivariate extreme-value copulas.

That part of the package has since been substantially redesigned.

The 1.0 extreme-value architecture is organized around the stable tail dependence function (STDF), with specialized bivariate Pickands-function machinery kept where it is useful.

That change allowed many models to move beyond dimension two, including Logistic, Galambos, Hüsler–Reiss, extremal-t, Tawn, asymmetric Galambos, Mixed, Cuadras–Augé, Marshall–Olkin, and discrete spectral constructions.

The same work generalized density and sampling infrastructure, introduced dimension-aware constructors, and added multivariate empirical extreme-value estimation.

So the extreme-value layer in 1.0 is not just “more families”; it is a different and more genuinely multivariate architecture.

Liouville copulas and real Williamson orders

Another addition that fits naturally with the generator/radial side of the package is Liouville copulas.

They support arbitrary positive Dirichlet parameters, arbitrary Archimedean generators under automatically checked monotonicity conditions, and integrate with the same high-level operations as the rest of the package: sampling, densities, CDFs, marginalization, conditioning, and Rosenblatt transforms.

Supporting Liouville models also pushed the Williamson machinery into more general territory, including non-integer orders and conditional radial representations.

This illustrates how new mathematical models have driven improvements in the common abstractions rather than being added as isolated implementations.

Student copulas are now much less of a special case

Student copulas have existed in Copulas.jl for a long time, but several numerical and inferential pieces were incomplete.

The work leading to 1.0 consolidated them:

• multivariate Student copula CDF evaluation;
• efficient scalar conditional distributions and h-functions;
• correct propagation of non-integer degrees of freedom under conditioning;
• bivariate Spearman’s rho;
• joint Kendall/Spearman rank matching;
• and shared scale-mixture probability machinery instead of separate numerical backends for related Student-based models.

A recurring theme of the 1.0 codebase is to share robust internal implementations whenever multiple models rely on the same mathematics.

Fitting became a real statistical-model interface

One of the largest changes since the earlier announcements is the fitting framework.

The original package already used Distributions.fit, but fitting has since grown into a much more coherent statistical-model API.

CopulaModel <: StatsBase.StatisticalModel represents a fitted copula or Sklar model. It retains the fitted distribution, the original fitting data, the fitted likelihood, and enough information to replay the estimator when a later procedure needs to refit the model. Coefficients, information criteria, null likelihoods, residuals, and related point-estimation summaries are exposed through the usual public statistical interfaces rather than through concrete fields.

One distinction became much sharper immediately before 1.0: fitting and post-fit uncertainty are separate operations.

Maximum likelihood (method=:mle) and maximum pseudo-likelihood (method=:mpl) are distinct procedures. Pseudo-observations have an explicit, order-invariant tie policy with configurable ranking methods. Rank-based estimators remain available where they make mathematical sense.

After fitting, uncertainty is computed explicitly with infer:

M = fit(CopulaModel, ClaytonCopula, U; method=:mle)
I = infer(M)

StatsBase.vcov(I)
StatsBase.stderror(I)
StatsBase.confint(I; level=0.95)

This keeps covariance information and resampling output out of the persistent CopulaModel state. Bootstrap and jackknife inference can instead replay the estimator recorded by the model, including sequential Sklar fits.

The same fitting interface also supports automatic family selection from an explicit candidate set:

Msel = fit(
    CopulaModel,
    Copulas.Copula,
    U;
    candidates=(ClaytonCopula, GumbelCopula, FrankCopula),
    criterion=:bic,
)

Mbest = selected_model(Msel)
selection_table(Msel)

AIC, BIC, AICc, and HQC are supported comparison criteria. Failed and non-finite fits remain visible in the selection table instead of silently disappearing.

The word “explicit” matters here. Copulas.jl does not pretend that every implemented family is a scientifically sensible candidate for every dataset and dimension. The user chooses the candidate set; the package provides a common fitting and comparison framework.

And the result objects themselves are deliberately opaque: retrieve the fitted distribution with fitted_distribution, the selected model with selected_model, and the comparison rows with selection_table. Their storage layout is not part of the 1.0 contract.

Parameter geometry without turning internals into API

The final stretch before 1.0 also involved a less visible refactor of parameter handling.

Copula families can have constrained scalar parameters, correlation matrices, nested generators, reflected models, Liouville parameters, or composite structures. Fitting all of these generically requires a notion of parameter geometry: how to validate a parameter point, move to unconstrained optimizer coordinates, and reconstruct a model afterwards.

That machinery is now delegated through Paramorph.jl rather than being reimplemented ad hoc across Copulas.jl.

This matters for maintainability and generic fitting, but it is intentionally not another user-facing protocol that downstream code must learn. The public contract stays semantic — constructors, params, fit, fitted_distribution, and the documented model operations — while optimizer coordinates and reconstruction details remain implementation machinery.

That separation is a useful example of what changed immediately before 1.0: the internals became more generic at the same time that the public surface became smaller and clearer.

From estimation to hypothesis testing

A statistical-model interface naturally raises the next question: how should those assumptions be tested?

Version 1.0 introduces a general copula hypothesis-testing framework built around a simple decomposition:

Hypothesis × Statistic × Calibration

The resulting CopulaTest machinery currently covers:

• independence;
• exchangeability;
• radial symmetry;
• extreme-value dependence;
• goodness of fit against a specified copula;
• goodness of fit for a fitted copula model.

Depending on the test, calibration can use simulation, randomization, multiplier bootstrap, or parametric bootstrap.

The design goal is the same as elsewhere in the package: adding a new hypothesis or a new statistic should primarily mean implementing the relevant mathematics, not inventing another unrelated user API.

Nonparametric dependence is first-class

Another part of the package that is easy to miss if one only remembers the early announcements is the nonparametric side.

Copulas.jl now includes EmpiricalCopula, BetaCopula, BernsteinCopula, CheckerboardCopula, empirical generators, and empirical extreme-value models.

They live in the same general distribution ecosystem as parametric copulas instead of being treated as a separate analysis module.

This matters for workflows where the goal is to estimate dependence without committing to a classical parametric family, or where a nonparametric fit is useful as a diagnostic/reference model.

Ecosystem interoperability kept expanding

The package began with Distributions.jl interoperability as its defining idea, and this has remained a useful architectural constraint.

Several optional integrations now build on top of it.

Plots.jl

A package extension provides recipes for copulas and SklarDist objects, including pairwise samples and bivariate CDF/density visualizations.

ExpectationMaximization.jl

Another extension supports weighted likelihood fitting and EM M-steps involving copulas, Sklar distributions, and compatible continuous mixture margins. This makes mixtures involving dependence models possible without changing the core fitting API.

PartitionedDistributions.jl

The newest integration connects the two packages’ marginalization and conditioning concepts.

For compatible vector-valued distributions that implement the required public marginal/conditional interface, Copulas.jl can expose subsetdims, condition, rosenblatt, and inverse_rosenblatt through that interface as well.

This is especially interesting architecturally because those operations are no longer confined to objects whose concrete type lives inside Copulas.jl.

Nataf correction

The package also provides Nataf, which computes the latent Gaussian correlation required for a Gaussian-copula/Sklar model to target a requested Pearson correlation matrix under specified non-Gaussian margins.

It is a small function compared with the larger fitting or testing frameworks, but a very practical one.

Probabilistic programming

And the old Turing story still holds: because Copula and SklarDist are distributions with likelihoods, they can participate naturally in probabilistic-programming models whenever the relevant numerical paths are differentiable.

That was one of the first compelling demonstrations of the design, and it remains one of the reasons to keep the package aligned with Julia’s standard statistical interfaces.

Numeric types are not an afterthought

The first announcement explicitly mentioned the goal of avoiding a hard dependency on Float64 throughout the package.

That principle has survived.

Not every algorithm can support every arbitrary real type, but across supported paths Copulas.jl tries to preserve Julia’s numeric genericity. Sampling and conditioning in particular now have explicit coverage for types such as Float32 and BigFloat.

There has also been a lot of less visible numerical work:

• generalized-quantile inversion that respects support boundaries, atoms, and plateaus;
• high-dimensional recurrences that avoid integer factorial overflow;
• stable inverse-derivative implementations for difficult Archimedean families;
• explicit boundary-density semantics;
• type-stable limiting constructors;
• lazy conditional computations;
• specialized dimension-aware fitting paths.

These changes are not headline features, but they are part of what makes a 1.0 release feel appropriate.

Performance: measure it, do not just claim it

Performance has always been one motivation for a native Julia implementation, but it is very easy to turn that into vague marketing.

Over time the package has added a proper Tachometer benchmark suite covering representative sampling, CDF/density, fitting, conditioning, and transform workloads. It is now used as an actual pre-release regression signal rather than only as a collection of performance examples. There is also a reproducible Julia/R comparison workflow.

The goal of these benchmarks is not to make the impossible claim that one language or package is universally faster. Different copula families and operations have very different computational structures.

The useful thing is that regressions are now measurable, and optimization work can be evaluated against concrete workloads.

Why 1.0 now?

A package can have a lot of features and still not be ready for 1.0 if users cannot tell which parts they may safely depend on. Conversely, 1.0 does not mean every imaginable copula family is implemented, every algorithm is optimal, or every future idea has already been designed. For Copulas.jl, the missing piece was an explicit public-API boundary.

The work immediately before 1.0 audited that boundary across the source code and documentation.

The rules are now intentionally simple:

1. Names declared export or public, together with their documented behavior and supported call forms, form the public Copulas.jl API.
2. Public documentation describes behavioral contracts, not accidental storage layouts.
3. Undocumented fields, caches, intermediate abstract hierarchies, dispatch traits, numerical backends, and optimization hooks remain internal.
4. The Manual, Bestiary, and Examples are narrative documentation; the Public API reference is the canonical behavioral reference.
5. Developer/internal documentation can explain implementation machinery without turning that machinery into a SemVer promise.

There are also a few deliberate extension contracts and exceptions. For example, Generator is intentionally extensible by downstream code. Its minimal mathematical interface is public.

Likewise, although arbitrary storage type parameters are not generally considered API, the form

SklarDist{CopulaType,Tuple{MarginTypes...}}

is deliberately supported as a fitting specification.

And for fitted models, users should retrieve the fitted distribution through fitted_distribution and use the public statistical-model and inference interfaces rather than depending on CopulaModel's concrete fields.

That level of precision is what makes the 1.0 designation appropriate.

Documentation as part of the API

The documentation has been reorganized around the same distinction.

The Manual should teach workflows and concepts.

The Bestiary should answer questions such as “what is this family?”, “what is its parameter domain?”, and “what dependence does it model?”.

The Examples should show realistic tasks.

The Public API should be the canonical reference for stable call forms, return semantics, ordering, dimensions, boundaries, and limitations.

The Developer Guide and Internal reference should help contributors understand how the implementation works without accidentally committing every internal hook to permanent compatibility.

This may look like documentation housekeeping, but for a 1.0 release it is part of the software design.

A short tour of Copulas.jl 1.0

A basic workflow remains intentionally ordinary:

using Copulas, Distributions

margins = (
    Gamma(2, 3),
    Beta(1, 4),
    Normal(),
)

C = ClaytonCopula(3, 2.0)
D = SklarDist(C, margins)

X = rand(D, 1_000)

From there, the same objects can participate in much richer workflows:

# Marginal/subset dependence model
C23 = subsetdims(C, (2, 3))

# Conditional distribution
Ccond = condition(C, 2, 0.4)

# Transform to independent uniforms and back
S = rosenblatt(D, X)
X2 = inverse_rosenblatt(D, S)

# Fit a full Sklar model
M = fit(
    CopulaModel,
    SklarDist{ClaytonCopula,Tuple{Gamma,Beta,Normal}},
    X,
)

Dfit = fitted_distribution(M)

And if the copula family is not known beforehand, the fitted-model framework can compare an explicit set of candidates using information criteria.

This composability is a more useful description of the package than simply counting the number of implemented copula families.

Compatibility and migration

Copulas.jl 1.0 requires:

Julia >= 1.11

If you are on an older Julia release, a compatible 0.1.x version will remain the appropriate choice.

For code migrating to 1.0, the main recommendation is to follow the now-explicit public interface:

• use documented constructors and methods;
• use public model accessors instead of concrete fields, and use infer for post-fit uncertainty;
• use fitted_distribution for fitted models;
• if you implement custom generators, rely on the documented Generator contract rather than internal derivative/fitting hooks;
• do not treat every concrete type parameter printed by Julia as a supported extension point.

The package has historically evolved quickly. The point of 1.0 is that future evolution should now be much less surprising for downstream users.

Looking back

The old Discourse thread also provides a useful record of how the project evolved.

In the first announcement, the project still lacked a proper testing suite and the roadmap listed things such as more Archimedean families, nested Archimedean copulas, Bernstein/Beta smoothing, and checkerboard constructions as future work.

Later, the community asked for nested copulas.

Then came the Williamson machinery, the JOSS review and publication, better documentation, bivariate extreme-value models, and many incremental releases whose changelogs gradually became longer than the public announcements.

Today the package has those nested models, the nonparametric families, broad fitting and post-fit inference, generalized conditioning and transforms, multivariate extreme-value models, Liouville copulas, hypothesis tests, model selection, ecosystem extensions, performance benchmarks, and a much more serious test architecture.

There are still plenty of things to improve and plenty of models that could be added. That is healthy.

The difference is that the project is no longer defined by whatever API the current implementation happens to expose.

It now has a contract.

Thank you

Copulas.jl 1.0 reflects contributions from many people.

Major pieces of the package have come from or been substantially improved by contributors including @Santymax98, @thisiscam, @FriesischScott, @langestefan, and many others through models, algorithms, tests, bug reports, documentation, and reviews. The project also benefited substantially from the JOSS review process, and the current logo was contributed by @lmiq.

Issues, difficult numerical examples, cross-checks against other implementations, feature requests, and API discussions have repeatedly exposed places where the package needed more general abstractions.

What comes after 1.0?

The 1.x series can continue in the same direction, but on a stable foundation.

The 1.x series can keep adding families, estimation procedures, statistical tests, numerical algorithms, performance improvements, and interoperability without requiring users to chase internal refactors.

That is the promise of this release.

For users working with dependence models in Julia, Copulas.jl 1.0 provides a stable base for sampling, estimation, inference, transformations, nonparametric models, and a broad range of copula families. Issues, discussions, benchmarks, mathematical references, and pull requests remain welcome.

Copulas.jl 1.0 is ready.

---

Links: Documentation · GitHub · Original Julia Discourse announcement thread