Add the rosenblatt transform and its inverse (#254)

Co-authored-by: Oskar Laverny <oskar.laverny@univ-amu.fr>

Author: FriesischScott

Project.toml

@@ -4,6 +4,7 @@ version = "0.1.30"
 authors = ["Oskar Laverny"]
 
 [deps]
+BigCombinatorics = "7b33fef7-ef1d-5d54-bb67-cf96d4c8a166"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -12,6 +13,7 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 MvNormalCDF = "37188c8d-bc69-4638-b057-733e744175ec"
+PolyLog = "85e3b03c-9856-11eb-0374-4dc1f8670e7f"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -23,6 +25,7 @@ WilliamsonTransforms = "48feb556-9bdd-43a2-8e10-96100ec25e22"
 
 [compat]
 Aqua = "v0.8"
+BigCombinatorics = "0.3.6"
 Combinatorics = "1"
 Distributions = "0.25"
 ForwardDiff = "0.10, 1"
@@ -32,6 +35,7 @@ InteractiveUtils = "1"
 LinearAlgebra = "1"
 LogExpFunctions = "0.3"
 MvNormalCDF = "0.2, 0.3"
+PolyLog = "2.6.0"
 PrecompileTools = "1"
 QuadGK = "2"
 Random = "1"
@@ -51,8 +55,9 @@ HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 TestItemRunner = "f8b46487-2199-4994-9208-9a1283c18c0a"
 
 [targets]
-test = ["Test", "TestItemRunner", "InteractiveUtils", "LinearAlgebra", "HypothesisTests", "Aqua", "StableRNGs"]
+test = ["Test", "TestItemRunner", "InteractiveUtils", "LinearAlgebra", "HypothesisTests", "Aqua", "StableRNGs", "StatsBase"]

docs/make.jl

@@ -31,6 +31,7 @@ makedocs(;
             "Empirical Copulas" => "empirical/generalities.md",
             "Vines Copulas" => "Vines.md",
             "Dependence measures" => "dependence_measures.md",
+            "Rosenblatt transformatons" =>  "rosenblatt_transform.md",
         ],
 
         "Bestiary" => [

docs/src/assets/references.bib

@@ -874,4 +874,15 @@ @article{blier2022stochastic
   pages={107506},
   year={2022},
   publisher={Elsevier}
+}
+
+@article{rosenblatt1952,
+  title={Remarks on a multivariate transformation},
+  author={Rosenblatt, Murray},
+  journal={The annals of mathematical statistics},
+  volume={23},
+  number={3},
+  pages={470--472},
+  year={1952},
+  publisher={JSTOR}
 }

docs/src/rosenblatt_transform.md

@@ -0,0 +1,45 @@
+```@meta
+CurrentModule = Copulas
+```
+
+# Rosenblatt transformations
+
+## Definition and usefullness
+
+
+!!! definition "Definition (Rosenblatt transformation):"
+    The Rosenblatt transformation considers a random vector ``X`` to be distributed as a certain multivariate cumulative distribution function ``F_{X}(x)``, and maps it back to a uniform distribution on the unit hypercube. 
+
+    More formally, consider the map ``R_X(x)`` defined as follows: 
+
+    ```math
+    R_X(x_1,...,x_d) = (r_1 = F_{X_1}(x_1), r_2 = F_{X_2 | X_1}(x_2 | x_1), ..., r_{d} = F_{X_d | X_1,...,X_{d-1}}(x_d|x_1,...x_{d-1}))
+    ```
+
+
+In certain circonstances, in paritcular for Archimedean copulas, this map simplifies to tractable expressions. it has a few nice properties: 
+
+* ``R_X(X) \sim \texttt{Uniform(Unit Hypercube)}``
+* ``R_X`` is a bijection. 
+
+These two properties are leveraged in some cases to construct the inverse rosenblatt transformations, that maps random noise to proper samples from the copula. In some cases, this is the best sampling algorithm available. 
+
+## Implementation
+
+As soon as the random vector ``X`` is represented by an object `X` that subtypes `SklarDist` or `Copula`, you have access to the `rosenblatt(X, x)` and `inverse_rosenblatt(X,x)` operators, which both have a straghtforward interpretation from their names. 
+
+```@docs
+rosenblatt
+inverse_rosenblatt
+```
+
+!!! note "Not all copulas available !"
+    Some copulas such has archimedeans have known expressions for their rosenblatt and/or inverse rosenblatt transforms, and therefore benefit from this interface and our implementation. On the other hand, some copulas have no known closed form expressions for conditional cdfs, and therefore their rosenblatt transformation is hard to implement.
+    
+    If you feel like you miss methods for certain particular copulas while the theory exists and it should be possible, do not hesitate to open an issue !
+
+
+```@bibliography
+Pages = [@__FILE__]
+Canonical = false
+```

src/ArchimedeanCopula.jl

@@ -1,33 +1,33 @@
 """
     ArchimedeanCopula{d, TG}
 
-Fields: 
-    - G::TG : the generator <: Generator. 
+Fields:
+    - G::TG : the generator <: Generator.
 
-Constructor: 
+Constructor:
 
     ArchimedeanCopula(d::Int,G::Generator)
 
-For some Archimedean [`Generator`](@ref) `G::Generator` and some dimenson `d`, this class models the archimedean copula which has this generator. The constructor checks for validity by ensuring that `max_monotony(G) ≥ d`. The ``d``-variate archimedean copula with generator ``\\phi`` writes: 
+For some Archimedean [`Generator`](@ref) `G::Generator` and some dimenson `d`, this class models the archimedean copula which has this generator. The constructor checks for validity by ensuring that `max_monotony(G) ≥ d`. The ``d``-variate archimedean copula with generator ``\\phi`` writes:
 
 ```math
 C(\\mathbf u) = \\phi^{-1}\\left(\\sum_{i=1}^d \\phi(u_i)\\right)
 ```
 
-The default sampling method is the Radial-simplex decomposition using the Williamson transformation of ``\\phi``. 
+The default sampling method is the Radial-simplex decomposition using the Williamson transformation of ``\\phi``.
 
-There exists several known parametric generators that are implement in the package. For every `NamedGenerator <: Generator` implemented in the package, we provide a type alias ``NamedCopula{d,...} = ArchimedeanCopula{d,NamedGenerator{...}}` to be able to manipulate the classic archimedean copulas without too much hassle for known and usefull special cases. 
+There exists several known parametric generators that are implement in the package. For every `NamedGenerator <: Generator` implemented in the package, we provide a type alias ``NamedCopula{d,...} = ArchimedeanCopula{d,NamedGenerator{...}}` to be able to manipulate the classic archimedean copulas without too much hassle for known and usefull special cases.
 
-A generic archimdean copula can be constructed as follows: 
+A generic archimdean copula can be constructed as follows:
 
 ```julia
 struct MyGenerator <: Generator end
 ϕ(G::MyGenerator,t) = exp(-t) # your archimedean generator, can be any d-monotonous function.
-max_monotony(G::MyGenerator) = Inf # could depend on generators parameters. 
+max_monotony(G::MyGenerator) = Inf # could depend on generators parameters.
 C = ArchimedeanCopula(d,MyGenerator())
 ```
 
-The obtained model can be used as follows: 
+The obtained model can be used as follows:
 ```julia
 spl = rand(C,1000)   # sampling
 cdf(C,spl)           # cdf
@@ -48,45 +48,34 @@ struct ArchimedeanCopula{d,TG} <: Copula{d}
         return new{d,typeof(G)}(G)
     end
 end
-ϕ(C::ArchimedeanCopula{d,TG},t)    where {d,TG} = ϕ(C.G,t)
-ϕ⁻¹(C::ArchimedeanCopula{d,TG},t)  where {d,TG} = ϕ⁻¹(C.G,t)
-ϕ⁽¹⁾(C::ArchimedeanCopula{d,TG},t) where {d,TG} = ϕ⁽¹⁾(C.G,t)
-ϕ⁽ᵏ⁾(C::ArchimedeanCopula{d,TG},t) where {d,TG} = ϕ⁽ᵏ⁾(C.G, Val(d), t)
-williamson_dist(C::ArchimedeanCopula{d,TG}) where {d,TG} = williamson_dist(C.G,Val(d))
+ϕ(C::ArchimedeanCopula{d,TG},t)      where {d,TG} = ϕ(C.G, t)
+ϕ⁻¹(C::ArchimedeanCopula{d,TG},t)    where {d,TG} = ϕ⁻¹(C.G, t)
+ϕ⁻¹⁽¹⁾(C::ArchimedeanCopula{d,TG},t) where {d,TG} = ϕ⁻¹⁽¹⁾(C.G, t)
+ϕ⁽¹⁾(C::ArchimedeanCopula{d,TG},t)   where {d,TG} = ϕ⁽¹⁾(C.G, t)
+ϕ⁽ᵏ⁾(C::ArchimedeanCopula{d,TG},k,t) where {d,TG} = ϕ⁽ᵏ⁾(C.G, k, t)
+williamson_dist(C::ArchimedeanCopula{d,TG}) where {d,TG} = williamson_dist(C.G, Val{d}())
 
-
-function _cdf(C::CT,u) where {CT<:ArchimedeanCopula} 
-    sum_ϕ⁻¹u = 0.0
-    for us in u
-        sum_ϕ⁻¹u += ϕ⁻¹(C,us)
-    end
-    return ϕ(C,sum_ϕ⁻¹u)
-end
+_cdf(C::ArchimedeanCopula, u) = ϕ(C, sum(ϕ⁻¹.(C, u)))
 function Distributions._logpdf(C::ArchimedeanCopula{d,TG}, u) where {d,TG}
-    if !all(0 .<= u .<= 1)
+    if !all(0 .< u .< 1)
         return eltype(u)(-Inf)
     end
-    T = promote_type(Float64, eltype(u)) # the Float64 here should be eltype(C) when copulas will be type agnostic... 
-    logdenom = sum_ϕ⁻¹u = zero(T)
-    for us in u
-        ϕ⁻¹u = ϕ⁻¹(C,us)
-        sum_ϕ⁻¹u += ϕ⁻¹u
-        logdenom += log(-ϕ⁽¹⁾(C,ϕ⁻¹u)) # log of negative here because ϕ⁽¹⁾ is necessarily negative
-    end
-    numer = abs(ϕ⁽ᵏ⁾(C, sum_ϕ⁻¹u))
-    if numer > 0
-        r = log(numer) - logdenom
-        !isnan(r) && return r
+    return log(ϕ⁽ᵏ⁾(C, Val{d}(), sum(ϕ⁻¹.(C, u))) * prod(ϕ⁻¹⁽¹⁾.(C, u)))
+end
+
+function Distributions._logpdf(C::ArchimedeanCopula{d,TG}, u) where {d,TG<:ClaytonGenerator}
+    if !all(0 .< u .< 1) || (C.G.θ < 0 && sum(u .^ -(C.G.θ)) < (d - 1))
+        return eltype(u)(-Inf)
     end
-    return -T(Inf)
+    return log(ϕ⁽ᵏ⁾(C, Val{d}(), sum(ϕ⁻¹.(C, u))) * prod(ϕ⁻¹⁽¹⁾.(C, u)))
 end
 
-# function τ(C::ArchimedeanCopula)  
+# function τ(C::ArchimedeanCopula)
 #     return 4*Distributions.expectation(r -> ϕ(C,r), williamson_dist(C)) - 1
 # end
 
 function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, x::AbstractVector{T}) where {T<:Real, CT<:ArchimedeanCopula}
-    # By default, we use the Williamson sampling. 
+    # By default, we use the Williamson sampling.
     Random.randexp!(rng,x)
     r = rand(rng,williamson_dist(C))
     sx = sum(x)
@@ -96,7 +85,7 @@ function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, x::Abstract
     return x
 end
 function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, A::DenseMatrix{T}) where {T<:Real, CT<:ArchimedeanCopula}
-    # More efficient version that precomputes the Williamson transform on each call to sample in batches: 
+    # More efficient version that precomputes the Williamson transform on each call to sample in batches:
     Random.randexp!(rng,A)
     n = size(A,2)
     r = rand(rng,williamson_dist(C),n)
@@ -112,7 +101,7 @@ function Distributions.fit(::Type{CT},u) where {CT <: ArchimedeanCopula}
     # @info "Archimedean fits are by default through inverse kendall tau."
     d = size(u,1)
     τ = StatsBase.corkendall(u')
-    # Then the off-diagonal elements of the matrix should be averaged: 
+    # Then the off-diagonal elements of the matrix should be averaged:
     avgτ = (sum(τ) .- d) / (d^2-d)
     GT = generatorof(CT)
     θ = τ⁻¹(GT,avgτ)
@@ -131,12 +120,12 @@ end
 ################################################################################################
 
 
-## Automatic syntactic sugar for all ZeroVariateGenerators and UnivariateGenerators. 
+## Automatic syntactic sugar for all ZeroVariateGenerators and UnivariateGenerators.
 ## see https://discourse.julialang.org/t/how-to-dispatch-on-a-type-alias/106476/38?u=lrnv
 function generatorof(::Type{S}) where {S <: ArchimedeanCopula}
     S2 = hasproperty(S,:body) ? S.body : S
     S3 = hasproperty(S2, :body) ? S2.body : S2
-    try 
+    try
         return S3.parameters[2].name.wrapper
     catch e
         @error "There is no generator type associated with the archimedean type $S"
@@ -159,7 +148,7 @@ for T in InteractiveUtils.subtypes(UnivariateGenerator)
     end
 end
 
-# The zero-variate ones just need a few more methods: 
+# The zero-variate ones just need a few more methods:
 Distributions._logpdf(::ArchimedeanCopula{d,IndependentGenerator}, u) where {d} = all(0 .<= u .<= 1) ? zero(eltype(u)) : eltype(u)(-Inf)
 Distributions._logpdf(::ArchimedeanCopula{d,MGenerator},           u) where {d} = all(u == u[1]) ?     zero(eltype(u)) : eltype(u)(-Inf)
 Distributions._logpdf(::ArchimedeanCopula{d,WGenerator},           u) where {d} = sum(u) == 1 ?   zero(eltype(u)) : eltype(u)(-Inf)
@@ -198,3 +187,76 @@ function Distributions._rand!(rng::Distributions.AbstractRNG, ::ArchimedeanCopul
     return A
 end
 
+function Distributions._rand!(rng::Distributions.AbstractRNG, C::ClaytonCopula, A::DenseMatrix{<:Real})
+    A[:] = inverse_rosenblatt(C, rand(rng, size(A)...))
+    return A
+end
+
+function rosenblatt(C::ArchimedeanCopula{d,TG}, u::AbstractMatrix{<:Real}) where {d,TG}
+    @assert d == size(u, 1)
+
+    U = zeros(eltype(u), size(u))
+    U[1, :] = u[1, :]
+    rⱼ = ϕ⁻¹.(C.G, u[1,:])
+    rⱼ₋₁ = similar(rⱼ)
+    for j in 2:d
+        rⱼ₋₁ .= rⱼ
+        rⱼ .+= ϕ⁻¹.(C.G, u[j,:]) # so we do not compute too much of them, nor allocate too much.
+        U[j, :] .= ϕ⁽ᵏ⁾.(C.G, Val(j - 1), rⱼ) ./ ϕ⁽ᵏ⁾.(C.G, Val(j - 1), rⱼ₋₁)
+    end
+    return U
+end
+
+function rosenblatt(
+    C::ArchimedeanCopula{d,TG}, u::AbstractMatrix{<:Real}
+) where {d,TG<:ClaytonGenerator}
+    @assert d == size(u, 1)
+
+    U = zeros(eltype(u), size(u))
+    U[1, :] = u[1, :]
+
+    for j in 2:d
+        U[j, :] .=
+            (
+                (1 .- j .+ sum(u[1:j, :] .^ (-C.G.θ); dims=1)[:]) ./
+                (2 .- j .+ sum(u[1:(j - 1), :] .^ (-C.G.θ); dims=1)[:])
+            ) .^ (-1 / C.G.θ - (j - 1))
+    end
+    return U
+end
+
+function inverse_rosenblatt(C::ArchimedeanCopula{d,TG}, u::AbstractMatrix{<:Real}) where {d,TG}
+    @assert d == size(u, 1)
+    U = zeros(eltype(u), size(u))
+    U[1, :] = u[1, :]
+    for i in axes(u, 2)
+        Cᵢⱼ = zero(eltype(u))
+        for j in 2:d
+            Cᵢⱼ += ϕ⁻¹(C.G, U[j - 1, i])
+            Dᵢⱼ = ϕ⁽ᵏ⁾(C.G, Val(j - 1), Cᵢⱼ) * u[j,i]
+            f(x) = ϕ⁽ᵏ⁾(C.G, Val(j - 1), Cᵢⱼ + ϕ⁻¹(C.G, x)) - Dᵢⱼ
+            U[j, i] = Roots.find_zero(f, (eps(1.0), 1.0), Roots.A42())
+        end
+    end
+    return U
+end
+
+function inverse_rosenblatt(
+    C::ArchimedeanCopula{d,TG}, u::AbstractMatrix{<:Real}
+) where {d,TG<:ClaytonGenerator}
+    @assert d == size(u, 1)
+
+    U = zeros(eltype(u), size(u))
+    U[1, :] = u[1, :]
+
+    for j in 2:d
+        U[j, :] =
+            (
+                1 .+
+                (1 .- (j - 1) .+ sum(U[1:(j - 1), :] .^ -C.G.θ; dims=1))[:] .*
+                (u[j, :] .^ (-1 / (j - 1 + 1 / C.G.θ)) .- 1)
+            ) .^ (-1 / C.G.θ)
+    end
+
+    return U
+end

src/BivariateArchimedeanMethods.jl

@@ -1,18 +1,18 @@
-# The Gumbel copula needs a specific bivariate method to handle large parameters. 
-function _cdf(C::ArchimedeanCopula{2,G}, u) where {G<:GumbelGenerator}
-    θ = C.G.θ
-    x₁, x₂ = -log(u[1]), -log(u[2])
-    lx₁, lx₂ = log(x₁), log(x₂)
-    return 1 - LogExpFunctions.cexpexp(LogExpFunctions.logaddexp(θ * lx₁, θ * lx₂) / θ)
-end
-function Distributions._logpdf(C::ArchimedeanCopula{2,G}, u) where {G<:GumbelGenerator}
-    T = promote_type(Float64, eltype(u))
-    !all(0 .< u .<= 1) && return T(-Inf) # if not in range return -Inf
-
-    θ = C.G.θ
-    x₁, x₂ = -log(u[1]), -log(u[2])
-    lx₁, lx₂ = log(x₁), log(x₂)
-    A = LogExpFunctions.logaddexp(θ * lx₁, θ * lx₂)
-    B = exp(A/θ)
-    return - B + x₁ + x₂ + (θ-1) * (lx₁ + lx₂) + A/θ - 2A + log(B + θ - 1)
-end
+# The Gumbel copula needs a specific bivariate method to handle large parameters.
+function _cdf(C::ArchimedeanCopula{2,G}, u) where {G<:GumbelGenerator}
+    θ = C.G.θ
+    x₁, x₂ = -log(u[1]), -log(u[2])
+    lx₁, lx₂ = log(x₁), log(x₂)
+    return 1 - LogExpFunctions.cexpexp(LogExpFunctions.logaddexp(θ * lx₁, θ * lx₂) / θ)
+end
+function Distributions._logpdf(C::ArchimedeanCopula{2,G}, u) where {G<:GumbelGenerator}
+    T = promote_type(Float64, eltype(u))
+    !all(0 .< u .<= 1) && return T(-Inf) # if not in range return -Inf
+
+    θ = C.G.θ
+    x₁, x₂ = -log(u[1]), -log(u[2])
+    lx₁, lx₂ = log(x₁), log(x₂)
+    A = LogExpFunctions.logaddexp(θ * lx₁, θ * lx₂)
+    B = exp(A/θ)
+    return - B + x₁ + x₂ + (θ-1) * (lx₁ + lx₂) + A/θ - 2A + log(B + θ - 1)
+end