# \[ANN\] Copulas.jl : A fully \`Distributions.jl\`-compliant copula package
**URL:** https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544
**Category:** Package Announcements
**Tags:** package, announcement, distributions, copula
**Created:** [February 16, 2022, 8:50am UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544 "2022-02-16T08:50:43Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [February 16, 2022, 8:50am UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/1 "2022-02-16T08:50:43Z")
**TLDR: I am very pleased to announce my first registered package: [`Copulas.jl`](https://github.com/lrnv/Copulas.jl) !**
> [@sylvaticus](#):
>
> For the unitialisated… [Copula (probability theory) - Wikipedia](https://en.wikipedia.org/wiki/Copula_(probability_theory))
> Also the links to the two leading packages in [R](https://cran.r-project.org/web/packages/copula/) and [Python](https://github.com/sdv-dev/Copulas)
## What is it
This package aims at bringing into native Julia most of the standard copula features: random number generation, fitting, construction of copula-based multivariate distributions through Sklar’s theorem, etc. while fully complying with the `Distributions.jl` API (after all, copulas are distributions functions).
Usually, people that use and work with copulas turn to R, because of the amazing `R` package `copula`.
While still perfectly maintained and updated today, the `R` package `copula` is full of obscured, heavily optimized, fast `C` code on one hand, and obscure, heavily optimized slow `R` code on the other hand (pun intended 😉 ).
This is an attempt to provide a very light, fast, reliable and maintainable copula implementation in native Julia (in particular, type-agnostic, so it’ll work with arbitrary type of floats like `Float32` for speed, `BigFloats` or `DoubleFloats` or `MultiFloats` for precision), with correct SIMD’sation, etc.
Two of the exported types are of most importance:
- `Copula` : this is an abstract mother type for all our copulas.
- `SklarDist` : Allows to construct a multivariate distribution by specifying the copula and the marginals, through Sklar’s theorem.
## What is already implemented
The API we implemented contains random number generations, cdf and pdf evaluations, and the `fit` function from `Distributions.jl`. Typical use case might look like this:
```julia
using Copulas, Distributions, Random
X₁ = Gamma(2,3)
X₂ = Pareto()
X₃ = LogNormal(0,1)
C = FranckCopula(3,7) # A 3-variate Franck Copula with θ = 7
D = SklarDist(C,(X₁,X₂,X₃)) # The final distribution
# This generates a (3,1000)-sized dataset from the multivariate distribution D
simu = rand(D,1000)
# While the following estimates the parameters of the model from a dataset :
D̂ = fit(SklarDist{ClaytonCopula,Tuple{Gamma,Normal,LogNormal}}, simu)
# Increase the number of observations to get a beter fit (well, might not work. can you spot why ?)
```
Atop from the very neat `SklarDist` type, available copulas are :
- `EmpiricalCopula`
- `GaussianCopula`
- `TCopula`
- `ArchimedeanCopula` (general, for any generator)
- `ClaytonCopula` (as an instantiation of `ArchimedeanCopula`, see after)
Next ones to be implemented will probably be :
- A few more archimedeans : `FranckCopula`, `AMHCopula`, `JoeCopula`, `GumbelCopula`
- Nested archimedeans (general, with the possibility to nest any family with any family, assuming it is possible, with parameter checks.)
- Bernstein copula and more general Beta copula as smoothing of the Empirical copula.
- `CheckerboardCopula` (and more generally `PatchworkCopula`)
Adding a new `ArchimedeanCopula` is very easy. The `Clayton` implementation is as short as :
```julia
struct ClaytonCopula{d,T} <: ArchimedeanCopula{d}
θ::T
end
ClaytonCopula(d,θ) = ClaytonCopula{d,typeof(θ)}(θ) # Constructor
radial_dist(C::ClaytonCopula) = Distributions.Gamma(1/C.θ,1) # Radial distribution
ϕ(C::ClaytonCopula, t) = (1+sign(C.θ)*t)^(-1/C.θ) # Generator
ϕ⁻¹(C::ClaytonCopula, t) = sign(C.θ)*(t^(-C.θ)-1) # Inverse Generator
τ(C::ClaytonCopula) = C.θ/(C.θ+2) # θ -> τ
τ⁻¹(::Type{ClaytonCopula},τ) = 2τ/(1-τ) # τ -> θ
```
The Archimedean API is modular :
- To sample an archimedean, only `radial_dist` and `ϕ` are needed.
- To evaluate the cdf, only `ϕ` and `ϕ⁻¹` are needed
- Currently to fit the copula, `τ⁻¹` is needed as we use inverted tau moment method. But we plan on also implementing inverse rho and MLE (density needed).
Both `FranckCopula`, `AMHCopula`, `JoeCopula`, `GumbelCopula` and many others are waiting for you to implement these methods !
## Dev roadmap
**With only a few dozens of lines of code, in a few hours, I reproduced a large part of the R `copula` package, which used thousands of lines of code for the same functionality…**
However, the release still lacks documentation and tests, which are on top of my todo-list. Nested archimedean copulas with very complient nesting ( **no more restriction to only nest from the same family** , as in R) is my next goal. Later, I’d like Bernstein & Beta copulas, Checkerboards & Patchworks, etc… The rest of the todo-list can be found [in the package’s README](https://github.com/lrnv/Copulas.jl#dev-roadmap).
I hope you’ll enjoy the package !
---
### Author:  [@rikh](https://discourse.julialang.org/u/rikh)
#### Post date: [February 16, 2022, 11:29am UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/2 "2022-02-16T11:29:07Z")
> [@lrnv](#):
>
> With only a few dozens of lines of code, in a few hours, I reproduced a large part of the R `copula` package, which used thousands of lines of code for the same functionality…
Exactly my experience with PowerAnalysis.jl versus R’s pwr package! I could avoid a lot of code duplication via Julia’s expressiveness. Julia is not all about raw throughput. Also about simplicity. Happy to hear that it worked out for Copulas.jl too
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [February 16, 2022, 11:42am UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/3 "2022-02-16T11:42:04Z")
> [@rikh](#):
>
> Exactly my experience with PowerAnalysis.jl versus R’s pwr package! I could avoid a lot of code duplication via Julia’s expressiveness. Julia is not all about raw passthrough. Also about simplicity. Happy to hear that it worked out for Copulas.jl too
The ease I had to code certain functionalities really surprised me, as I remember coding some `S4` classes in R (in another life) that were really cumbersome… When I started I thought this was a project for the summer, and in a week-end I did half of it. I realize now that using the `Distributions.jl` framework is actually what allowed me to do this as efficiently.
---
### Author:  [@juliohm](https://discourse.julialang.org/u/juliohm)
#### Post date: [February 16, 2022, 11:43am UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/4 "2022-02-16T11:43:37Z")
Super nice to know there is a Copulas.jl package now. Thanks for the contribution @lrnv!
> [@rikh](#):
>
> Exactly my experience with PowerAnalysis.jl versus R’s pwr package!
I had the exact same experience writing CoDa.jl, the equivalent of R’s compositions package. It is so nice to code in Julia!
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [February 16, 2022, 11:45am UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/5 "2022-02-16T11:45:58Z")
> [@juliohm](#):
>
> Super nice to know there is a Copulas.jl package now. Thanks for the contribution @lrnv!
Thanks ! As i said, this is only a first draft and many things are lacking in the package… starting with a propper testing suite 😉 But my goal is parity with the R package and even beyond.
---
### Author:  [@tbeason](https://discourse.julialang.org/u/tbeason)
#### Post date: [February 16, 2022, 1:49pm UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/6 "2022-02-16T13:49:36Z")
This is great! I’ll be sure to check it out in the near future.
---
### Author:  [@sylvaticus](https://discourse.julialang.org/u/sylvaticus)
#### Post date: [February 16, 2022, 1:57pm UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/7 "2022-02-16T13:57:25Z")
For the unitialisated… [Copula (probability theory) - Wikipedia](https://en.wikipedia.org/wiki/Copula_(probability_theory))
Also the links to the two leading packages in [R](https://cran.r-project.org/web/packages/copula/) and [Python](https://github.com/sdv-dev/Copulas)
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [February 16, 2022, 2:20pm UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/8 "2022-02-16T14:20:38Z")
@sylvaticus quoting you in Op, hope you dont mind.
---
### Author:  [@xiaodai](https://discourse.julialang.org/u/xiaodai)
#### Post date: [February 17, 2022, 1:21am UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/9 "2022-02-17T01:21:38Z")
nice one
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [February 21, 2022, 11:07am UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/10 "2022-02-21T11:07:42Z")
#### Version 0.1.1 is out
This first patch version brings a few more copulas to the table : `FranckCopula` , `AMHCopula` , `JoeCopula` , `GumbelCopula`. All of them use the Archimedean Frailty representation for simulations.
The documentation now contains the definition of all exported copulas.
Next step : testing 🙂
---
### Author:  [@AAL](https://discourse.julialang.org/u/AAL)
#### Post date: [February 22, 2022, 2:07pm UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/11 "2022-02-22T14:07:48Z")
I use the [DatagenCopulaBased](https://github.com/iitis/DatagenCopulaBased.jl) package in some of my projects. @lrnv You seem to be also actively involved in that project. Do you recommend switching to this package? Are there any major differences between the two packages I should be aware of before switching?
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [February 22, 2022, 2:24pm UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/12 "2022-02-22T14:24:23Z")
> [@AAL](#):
>
> I use the [DatagenCopulaBased](https://github.com/iitis/DatagenCopulaBased.jl) package in some of my projects. @lrnv You seem to be also actively involved in that project. Do you recommend switching to this package? Are there any major differences between the two packages I should be aware of before switching?
I’m not much involved, I just fixed a few things that were bothering me the other week and then decided that a more general package might be needed. However:
- `DatagenCopulaBased.jl` is for the moment much better tested than `Copulas.jl`, therefore I’ll recommend to keep using it for now (and open issues / PR onto `Copulas.jl` for missing features if you have time / want, you’ll always be welcome 🙂 ). Hopefully in the few next months, the testing and feature levels of `Copulas.jl` will increase and allow you to switch over.
- The main difference between the two is that `DatagenCopulaBased.jl`, as its name suggests, only provides random number generations, while `Copulas.jl` provides full interactivity with `Distributions.jl`’s methods, which allows to use `Distributions.fit()`, and possibility to interact with the broader ecosystem such as `Turing.jl`, `StatsPlots.jl` etc. There is also the new `SklarDist{}` type, which is also a `Distributions.jl`-complient type which is AFAIK completely new in the ecosystem.
* * *
PS: There are also more copulas available in `DatagenCopulaBased`. One of my goals is indeed parity, which will come someday, some of the code will be ported from one package to the other soon.
* * *
PPS: The `SklarDist` construct allows you to do stuff as :
```julia
julia> MyDist = SklarDist(SurvivalCopula(FrankCopula(4,7),(2,4)),(LogNormal(),Pareto(),Gamma(),Normal()));
julia> data = rand(MyDist,1000)
4×1000 Matrix{Float64}:
2.35294 0.457584 2.44058 0.840622 1.81306 0.630984 … 2.22521 0.496158 0.525926 0.701049 0.87533 4.92817
1.47767 42.4113 1.55615 14.7199 1.75489 1.01644 13.6979 80.3027 1.40588 1.51737 1.40788 1.16863
3.40522 0.0811707 2.19233 0.154554 2.5537 0.342503 0.0466783 0.134426 0.943388 0.853033 0.557026 1.49137
-0.0871605 1.1993 -2.18685 0.626663 -0.553239 0.133756 0.668708 0.81423 0.271098 -0.295031 0.00926165 -0.891382
julia> MyCop = SurvivalCopula{4,ClaytonCopula,(2,4)};
julia> MyMarginals = Tuple{LogNormal,Pareto,Gamma,Normal};
julia> fitted_model = fit(SklarDist{MyCop,MyMarginals},data)
```
Which works correctly, and estimates the Copula, taking into acount the fact that it is `Survival` in dimensions 2 and 4, (meaning that the pseudo-observations are flipped), together with the margins. I think `SklarDist{Cop,Margins}` is the right name for this distribution.
You could also prefer to keep the empirical copula in the multivariate distribution, which you can easily do via :
```julia
julia> MyCop = EmpiricalCopula;
julia> fitted_model = fit(SklarDist{MyCop,MyMarginals},data)
```
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [June 20, 2022, 11:02am UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/13 "2022-06-20T11:02:29Z")
**Version 0.1.4 is out.**
It contains only a few bug-fixes, and one new feature : implementation of a generic algorithm for Archimedean densities. This allows to leverage `Distributions.jl`’s compatability with the broader ecosystem, i.e. to do MCMC with `Turing.jl` on complicated models as follows:
```julia
using Copulas, Distributions, Random, Turing, Plots, StatsPlots
Random.seed!(123)
D = SklarDist(ClaytonCopula(3,7), (Exponential(1.0), Pareto(3.0), Exponential(2.0)))
draws = rand(D, 2_000)
@model function copula(X)
# Priors
θ ~ TruncatedNormal(1.0, 1.0, 0, Inf)
θ₁ ~ TruncatedNormal(1.0, 1.0, 0, Inf)
θ₂ ~ TruncatedNormal(1.0, 1.0, 0, Inf)
θ₃ ~ TruncatedNormal(1.0, 1.0, 0, Inf)
# Marginals distributions and copula
X₁ = Exponential(θ₁)
X₂ = Pareto(θ₂)
X₃ = Exponential(θ₃)
C = ClaytonCopula(3,θ)
D = SklarDist(C, (X₁, X₂, X₃))
Turing.Turing.@addlogprob! loglikelihood(D, X)
end
sampler = NUTS() # MH() works too
chain = sample(copula(draws), sampler, MCMCThreads(), 1_00, 4)
p = plot(chain)
savefig(p, "plot_archimedeans.png")
```
---
### Author:  [@Domenic\_Di\_Francesco](https://discourse.julialang.org/u/Domenic_Di_Francesco)
#### Post date: [October 17, 2022, 12:48pm UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/14 "2022-10-17T12:48:21Z")
@lrnv - the Turing compatibility in particular looks great!
How can I download version 0.1.4 of your package? The Julia package manager suggests that v0.1.0 is the latest available.
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [October 17, 2022, 1:50pm UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/15 "2022-10-17T13:50:27Z")
@Domenic_Di_Francesco Thanks, glad you enjoy it !
This is weird… Did you try `] up` ? I think it might be a compatibility issue in your environment.
The most up to date version is 0.1.6, and the Julia package manager should be able to install it. Updating your packages should be enough to resolve the compatibility issue
---
### Author:  [@Domenic\_Di\_Francesco](https://discourse.julialang.org/u/Domenic_Di_Francesco)
#### Post date: [October 17, 2022, 2:53pm UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/16 "2022-10-17T14:53:46Z")
@lrnv - hi, just to confirm that after various updates the issue is resolved for me.
Thank you again for this contribution! The straightforward integration with `Distributions` and `Turing` is fantastic.
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [October 17, 2022, 3:01pm UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/17 "2022-10-17T15:01:07Z")
> [@Domenic\_Di\_Francesco](#):
>
> @lrnv - hi, just to confirm that after various updates the issue is resolved for me.
Glad to hear it then.
> [@Domenic\_Di\_Francesco](#):
>
> Thank you again for this contribution! The straightforward integration with `Distributions` and `Turing` is fantastic.
Indeed ! The best part is the ridiculous amount of code that was needed for that to happen.
---
### Author:  [@Kokora](https://discourse.julialang.org/u/Kokora)
#### Post date: [November 26, 2022, 9:52pm UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/18 "2022-11-26T21:52:18Z")
@lrnv thank you for your contribution on this. This is great!
Is it possible to implement a regression version of this? Namely, how can I specify a regression model for the θ, θ₁, θ₂?
I created a small example (see below). But it seems to not work very well. I am not sure what I am doing wrong. Any help would be greatly appreciated.
```julia
using Random
Random.seed!(43);
n = 200
p = 2
X = randn(n, p)
β1 = randn(p)
Y = X*hcat(β1,β1) + transpose(rand(SklarDist(GaussianCopula([1 .8; .8 1]), (Normal(0.0, 1.), Normal(0.0,1.0))), n))#transpose(rand(MultivariateNormal(zeros(2), [1 0.8; 0.8 1]),n))
#-- LOok at the data
cornerplot(Y, compact=true)
#-- Model
@model function copulaReg(Y, X)
n, K = size(X)
# Priors
ρ ~ truncated(Normal(0,.3), -1, 1) #-- Corelated
θ ~ Normal(0.,1.)
θ₁ ~ Normal(0.,1.)
θ₂ ~ Normal(0.,1.)
β1 ~ filldist(Normal(0., 1.), K)
β2 ~ filldist(Normal(0., 1.), K)
β ~ filldist(Normal(0., 1.), K)
#-- Regression parameters
μ1 = zeros(n)
μ2 = zeros(n)
for i in 1:n
μ1[i] = θ₁ .+ sum(X[i,:].* β1)
μ2[i] = θ₂ .+ sum(X[i,:].* β2)
#Turing.Turing.@addlogprob! loglikelihood(SklarDist(GaussianCopula([1. ρ; ρ 1.0]), (Normal(μ1[i], σ1), Normal(μ2[i], σ2))), Y[i,:])
Turing.Turing.@addlogprob! loglikelihood(SklarDist(GaussianCopula([1. ρ; ρ 1.0]), (Normal(μ1[i], 1.), Normal(μ2[i], 1.))), Y[i,:])
end
end
#-- sampler
sampler = NUTS() # doe not work
sampler = PG(10) # works but gives weird results!!
#-- sample from the posterior
chain = sample(copulaReg(Y, X), sampler, MCMCThreads(), 2_000, 3)
describe(chain)
```
---
### Author:  [@lrnv](https://discourse.julialang.org/u/lrnv)
#### Post date: [November 27, 2022, 2:48pm UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/19 "2022-11-27T14:48:18Z")
@Kokora Yes, this regression version should work out of the box.
Unfortunately, I tried everything I could and still gets the same error every time. Could you please open a topic specific to this problem in the [Specific Domains / Probabilistic programming](https://discourse.julialang.org/c/domain/probprog/48) so that we could draw attention from people that might know more than me the internals of `Turing.jl` ? That would be great.
---
### Author:  [@Kokora](https://discourse.julialang.org/u/Kokora)
#### Post date: [November 27, 2022, 3:30pm UTC](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544/20 "2022-11-27T15:30:32Z")
@lrn Thank you for looking at this. I will go ahead create a topic specific as you suggested.
[Next page](https://discourse.julialang.org/t/ann-copulas-jl-a-fully-distributions-jl-compliant-copula-package/76544.md?page=2)