Wiliamson transform sampling v2

Project.toml

@@ -14,19 +14,21 @@ Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
+WilliamsonTransforms = "48feb556-9bdd-43a2-8e10-96100ec25e22"
 
 [compat]
 Cubature = "1.5"
 Distributions = "0.25"
 ForwardDiff = "0.10"
 GSL = "1"
 MvNormalCDF = "0.2, 0.3"
+Random = "1.6"
 Roots = "1, 2"
 SpecialFunctions = "2"
 StatsBase = "0.33, 0.34"
 TaylorSeries = "0.12, 0.13, 0.14, 0.15"
+WilliamsonTransforms = "0.1"
 julia = "1.6"
-Random = "1.6"
 
 [extras]
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

README.md

@@ -41,16 +41,17 @@ simu = rand(D,1000)
 D̂ = fit(SklarDist{FrankCopula,Tuple{Gamma,Normal,LogNormal}}, simu)
 # Increase the number of observations to get a beter fit (or not?)  
 ```
+
 Available copula families are:
-- `GaussianCopula`,
-- `TCopula`,
-- `ArchimedeanCopula` (for any generator),
-- `ClaytonCopula`,`FrankCopula`, `AMHCopula`, `JoeCopula`, `GumbelCopula` as example of the `ArchimedeanCopula` abstract type, see below,
-- `WCopula` and `MCopula`, which are [Fréchet-Hoeffding bounds](https://en.wikipedia.org/wiki/Copula_(probability_theory)#Fr%C3%A9chet%E2%80%93Hoeffding_copula_bounds),
+- `EllipticalCopulas`: `GaussianCopula` and `TCopula`
+- `ArchimedeanCopula`: `WilliamsonCopula` (for any generator), but also `ClaytonCopula`,`FrankCopula`, `AMHCopula`, `JoeCopula`, `GumbelCopula`, supporting the full ranges in every dimensions (e.g. ClaytonCopula can be sampled with negative dependence in any dimension, not just d=2). 
+- `WCopula`, `IndependentCopula` and `MCopula`, which are [Fréchet-Hoeffding bounds](https://en.wikipedia.org/wiki/Copula_(probability_theory)#Fr%C3%A9chet%E2%80%93Hoeffding_copula_bounds),
+- `PlackettCopula`, see ref?
 - `EmpiricalCopula` to follow closely a given dataset.
 
 The next ones to be implemented will probably be: 
-- Nested archimedeans (general, with the possibility to nest any family with any family, assuming it is possible, with parameter checks.)
+- Extreme values copulas. 
+- Nested archimedeans (for any generators, with automatic nesting conditions checking). 
 - Bernstein copula and more general Beta copula as smoothing of the Empirical copula. 
 - `CheckerboardCopula` (and more generally `PatchworkCopula`)
 
@@ -65,35 +66,30 @@ ClaytonCopula(d,θ)            = ClaytonCopula{d,typeof(θ)}(θ)     # Construct
 ϕ⁻¹(C::ClaytonCopula,t)       = sign(C.θ)*(t^(-C.θ)-1)            # Inverse Generator
 τ(C::ClaytonCopula)           = C.θ/(C.θ+2)                       # θ -> τ
 τ⁻¹(::Type{ClaytonCopula},τ)  = 2τ/(1-τ)                          # τ -> θ
-radial_dist(C::ClaytonCopula) = Distributions.Gamma(1/C.θ,1)      # Radial distribution
+williamson_dist(C::ClaytonCopula{d,T}) where {d,T} = WilliamsonFromFrailty(Distributions.Gamma(1/C.θ,1),d) # Radial distribution
 ```
 The Archimedean API is modular: 
 
-- To sample an archimedean, only `radial_dist` and `ϕ` are needed.
+- To sample an archimedean, only `ϕ` is required. Indeed, the `williamson_dist` has a generic fallback that uses [WilliamsonTransforms.jl](https://www.github.com/lrnv/WilliamsonTransforms.jl) for any generator. Note however that providing the `williamson_dist` yourself if you know it will allow sampling to be an order of magnitude faster.
 - To evaluate the cdf and (log-)density in any dimension, only `ϕ` and `ϕ⁻¹` are needed.
 - Currently, to fit the copula `τ⁻¹` is needed as we use the inverse tau moment method. But we plan on also implementing inverse rho and MLE (density needed). 
 - Note that the generator `ϕ` follows the convention `ϕ(0)=1`, while others (e.g., https://en.wikipedia.org/wiki/Copula_(probability_theory)#Archimedean_copulas) use `ϕ⁻¹` as the generator.
-- We plan on implementing the Williamson transformations so that `radial-dist` can be automaticlaly deduced from `ϕ` and vice versa, if you dont know much about your archimedean family
 
 # Dev Roadmap
 
-## Urgent things
-
-- [ ] Add tests and documentation
-
 ## Next
 
-- [ ] Extensive documentation and tests for the current implementation. 
-- [x] Implement archimedean density generally. 
-- [ ] Docs: show how to implement another archimedean.  
+- [ ] More documentation and tests for the current implementation. 
+- [ ] Docs: show how to use the WilliamsonCopula to implement generic archimedeans.
 - [ ] Give the user the choice of fitting method via `fit(dist,data; method="MLE")` or `fit(dist,data; method="itau")` or `fit(dist,data; method="irho")`.
-- [ ] Fitting a generic archimedean : should provide an empirical generator
-- [ ] Make `Archimedean` more generic : inputing only `radial_dist` or only `phi` shoudl be enough to get `pdf, cdf, rand, tau, rho, itau, irho, fit, radial_dist`, etc...  **Williamson d-transform and inverse d-transform should be implemented.** The checking of nesting possibility should be done automatically with some rules (is phi_inv \circ phi complementely monotonous ? with obviously shortcut for inter-family nestings.)   
+- [ ] Fitting a generic archimedean with an empirically produced generarator
+- [ ] Automatic checking of generator d-monotonicity ? Dunno if it is even possible. 
 
 ## Maybe later
 
+- [ ] `NestedArchimedean`, with automatic checking of nesting conditions for generators. 
 - [ ] `Vines`?
-- [ ] `NestedArchimedean` and very easy implementation of new archimeean copulas via the radial dist or the phi/invphi + Williamson transform. 
+- [ ] `Archimax` ?
 - [ ] `BernsteinCopula` and `BetaCopula` could also be implemented. 
 - [ ] `PatchworkCopula` and `CheckerboardCopula`: could be nice things to have :)
 - [ ] Goodness of fits tests ?

docs/src/archimedean/generalities.md

@@ -4,10 +4,28 @@ CurrentModule = Copulas
 
 # General discussion
 
-Details about what are archimedean copulas (on the math level)
+Archimedean copulas are a large class of copulas that are defined as : ... [ details needed...]
+
+Adding a new `ArchimedeanCopula` is very easy. The `Clayton` implementation is as short as: 
+
+```julia
+struct ClaytonCopula{d,T} <: Copulas.ArchimedeanCopula{d}
+    θ::T
+end
+ClaytonCopula(d,θ)            = ClaytonCopula{d,typeof(θ)}(θ)     # Constructor
+ϕ(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-τ)                          # τ -> θ
+williamson_dist(C::ClaytonCopula{d,T}) where {d,T} = WilliamsonFromFrailty(Distributions.Gamma(1/C.θ,1),d) # Radial distribution
+```
+The Archimedean API is modular: 
+
+- To sample an archimedean, only `ϕ` is required. Indeed, the `williamson_dist` has a generic fallback that uses [WilliamsonTransforms.jl](https://www.github.com/lrnv/WilliamsonTransforms.jl) for any generator. Note however that providing the `williamson_dist` yourself if you know it will allow sampling to be an order of magnitude faster.
+- To evaluate the cdf and (log-)density in any dimension, only `ϕ` and `ϕ⁻¹` are needed.
+- Currently, to fit the copula `τ⁻¹` is needed as we use the inverse tau moment method. But we plan on also implementing inverse rho and MLE (density needed). 
+- Note that the generator `ϕ` follows the convention `ϕ(0)=1`, while others (e.g., https://en.wikipedia.org/wiki/Copula_(probability_theory)#Archimedean_copulas) use `ϕ⁻¹` as the generator.
 
-Details about the generic construction of archimedean copulas (in the package), 
-and details on exported methods that corresponds to this class 
 
 ```@docs
 ArchimedeanCopula

docs/src/archimedean/implement_your_own.md

@@ -14,5 +14,5 @@ end
 ϕ⁻¹(C::MyAC{d},x)
 τ(C::MyAC{d})
 τ⁻¹(::MyAC{d},τ)
-radial_dist(C::MyAC{d})
+williamson_dist(C::MyAC{d})
 ```

src/ArchimedeanCopula.jl

@@ -26,15 +26,18 @@ ClaytonCopula(d,θ)            = ClaytonCopula{d,typeof(θ)}(θ)     # Construct
 ϕ⁻¹(C::ClaytonCopula,t)       = sign(C.θ)*(t^(-C.θ)-1)            # Inverse Generator
 τ(C::ClaytonCopula)           = C.θ/(C.θ+2)                       # θ -> τ
 τ⁻¹(::Type{ClaytonCopula},τ)  = 2τ/(1-τ)                          # τ -> θ
-radial_dist(C::ClaytonCopula) = Distributions.Gamma(1/C.θ,1)      # Radial distribution
+williamson_dist(C::ClaytonCopula{d,T}) where {d,T} = WilliamsonFromFrailty(Distributions.Gamma(1/C.θ,1),d) # Radial distribution
 ```
 The Archimedean API is modular: 
 
-- To sample an archimedean, only `radial_dist` and `ϕ` are needed.
+- To sample an archimedean, only `williamson_dist` and `ϕ` are needed.
 - To evaluate the cdf and (log-)density in any dimension, only `ϕ` and `ϕ⁻¹` are needed.
 - Currently, to fit the copula `τ⁻¹` is needed as we use the inverse tau moment method. But we plan on also implementing inverse rho and MLE (density needed). 
 - Note that the generator `ϕ` follows the convention `ϕ(0)=1`, while others (e.g., https://en.wikipedia.org/wiki/Copula_(probability_theory)#Archimedean_copulas) use `ϕ⁻¹` as the generator.
 - We plan on implementing the Williamson transformations so that `radial-dist` can be automaticlaly deduced from `ϕ` and vice versa, if you dont know much about your archimedean family
+
+If you only know the generator of your copula, take a look at WilliamsonCopula that allows to generate automatically the associated williamson distribution. 
+If on the other hand you have a univaraite positive random variable with no atom at zero, then the williamson transform can produce an archimdean copula out of it, with the same constructor. 
 """
 abstract type ArchimedeanCopula{d} <: Copula{d} end
 function Distributions.cdf(C::CT,u) where {CT<:ArchimedeanCopula} 
@@ -64,26 +67,42 @@ function Distributions._logpdf(C::CT, u) where {CT<:ArchimedeanCopula}
     for us in u
         ϕ⁻¹u = ϕ⁻¹(C,us)
         sum_ϕ⁻¹u += ϕ⁻¹u
-        sum_logϕ⁽¹⁾ϕ⁻¹u += log(-ϕ⁽¹⁾(C,ϕ⁻¹u)) #log of negative here because ϕ⁽¹⁾ is necessarily negative
+        sum_logϕ⁽¹⁾ϕ⁻¹u += log(-ϕ⁽¹⁾(C,ϕ⁻¹u)) # log of negative here because ϕ⁽¹⁾ is necessarily negative
     end
-
     numer = ϕ⁽ᵈ⁾(C, sum_ϕ⁻¹u)
+    @show sum_logϕ⁽¹⁾ϕ⁻¹u, sum_ϕ⁻¹u, numer
     dimension_sign = iseven(d) ? 1.0 : -1.0 #need this for log since (-1.0)ᵈ ϕ⁽ᵈ⁾ ≥ 0.0
-    return log(dimension_sign*numer) - sum_logϕ⁽¹⁾ϕ⁻¹u
+
+
+    # I am not sure this is the right reasoning :
+    if numer == 0
+        if sum_logϕ⁽¹⁾ϕ⁻¹u == -Inf
+            return Inf
+        else
+            return -Inf
+        end
+    else
+        return log(dimension_sign*numer) - sum_logϕ⁽¹⁾ϕ⁻¹u
+    end
 end
 
 ϕ(C::ArchimedeanCopula{d},x) where d = @error "Archimedean interface not implemented for $(typeof(C)) yet."
 ϕ⁻¹(C::ArchimedeanCopula{d},x) where d = @error "Archimedean interface not implemented for $(typeof(C)) yet."
-radial_dist(C::ArchimedeanCopula{d}) where d = @error "Archimedean interface not implemented for $(typeof(C)) yet."
 τ(C::ArchimedeanCopula{d}) where d  = @error "Archimedean interface not implemented for $(typeof(C)) yet."
 τ⁻¹(::ArchimedeanCopula{d},τ) where d = @error "Archimedean interface not implemented for $(typeof(C)) yet."
-# radial_dist(C::ArchimedeanCopula) = laplace_transform(t -> ϕ(C,t))
+williamson_dist(C::ArchimedeanCopula{d}) where d = WilliamsonTransforms.𝒲₋₁(t -> ϕ(C,t),d)
 
 function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, x::AbstractVector{T}) where {T<:Real, CT<:ArchimedeanCopula}
-    r = rand(rng,radial_dist(C))
+    # By default, we use the williamson sampling. 
+    # the copula must have a field R that correspond to its williamson transformed variable for this to work.
     Random.rand!(rng,x)
+    r = rand(rng,williamson_dist(C))
     for i in 1:length(C)
-        x[i] = ϕ(C,-log(x[i])/r)
+        x[i] = -log(x[i])
+    end
+    sx = sum(x)
+    for i in 1:length(C)
+        x[i] = ϕ(C,r * x[i]/sx)
     end
     return x
 end
@@ -95,4 +114,4 @@ function Distributions.fit(::Type{CT},u) where {CT <: ArchimedeanCopula}
     avgτ = (sum(τ) .- d) / (d^2-d)
     θ = τ⁻¹(CT,avgτ)
     return CT(d,θ)
-end
+end

src/ArchimedeanCopulas/AMHCopula.jl

@@ -21,7 +21,7 @@ It has a few special cases:
 struct AMHCopula{d,T} <: ArchimedeanCopula{d}
     θ::T
     function AMHCopula(d,θ)
-        if (θ < -1) || (θ > 1)
+        if (θ < -1) || (θ >= 1)
             throw(ArgumentError("Theta must be in [-1,1)"))
         elseif θ == 0
             return IndependentCopula(d)
@@ -44,8 +44,6 @@ function τ⁻¹(::Type{AMHCopula},τ)
     end
     return Roots.fzero(θ -> 1 - 2(θ+(1-θ)^2*log(1-θ))/(3θ^2) - τ,0.5)
 end
-
-
-radial_dist(C::AMHCopula) = 1 + Distributions.Geometric(1-C.θ)
+williamson_dist(C::AMHCopula{d,T}) where {d,T} = C.θ >= 0 ? WilliamsonFromFrailty(1 + Distributions.Geometric(1-C.θ),d) : WilliamsonTransforms.𝒲₋₁(t -> ϕ(C,t),d)
 
 

src/ArchimedeanCopulas/ClaytonCopula.jl

@@ -36,8 +36,8 @@ struct ClaytonCopula{d,T} <: ArchimedeanCopula{d}
         end
     end
 end
-ϕ(  C::ClaytonCopula,      t) = (1+sign(C.θ)*t)^(-1/C.θ)
-ϕ⁻¹(C::ClaytonCopula,      t) = sign(C.θ)*(t^(-C.θ)-1)
-radial_dist(C::ClaytonCopula) = Distributions.Gamma(1/C.θ,1) # Currently fails for negative thetas ! thus negtatively correlated clayton copulas cannot be sampled...
+ϕ(  C::ClaytonCopula,      t) = max(1+C.θ*t,zero(t))^(-1/C.θ)
+ϕ⁻¹(C::ClaytonCopula,      t) = (t^(-C.θ)-1)/C.θ
 τ(C::ClaytonCopula) = ifelse(isfinite(C.θ), C.θ/(C.θ+2), 1)
-τ⁻¹(::Type{ClaytonCopula},τ) = ifelse(τ == 1,Inf,2τ/(1-τ))
+τ⁻¹(::Type{ClaytonCopula},τ) = ifelse(τ == 1,Inf,2τ/(1-τ))
+williamson_dist(C::ClaytonCopula{d,T}) where {d,T} = C.θ >= 0 ? WilliamsonFromFrailty(Distributions.Gamma(1/C.θ,1),d) : ClaytonWilliamsonDistribution(C.θ,d)

src/ArchimedeanCopulas/FrankCopula.jl

@@ -55,7 +55,6 @@ function τ⁻¹(::Type{FrankCopula},τ)
     return Roots.fzero(x -> (1-D₁(x))/x - x₀, 1e-4, Inf)
 end
 
-
-radial_dist(C::FrankCopula) = Logarithmic(1-exp(-C.θ))
+williamson_dist(C::FrankCopula{d,T}) where {d,T} = WilliamsonFromFrailty(Logarithmic(1-exp(-C.θ)), d)
 
 

src/ArchimedeanCopulas/GumbelCopula.jl

@@ -37,8 +37,7 @@ end
 ϕ⁻¹(C::GumbelCopula,       t) = (-log(t))^C.θ
 τ(C::GumbelCopula) = ifelse(isfinite(C.θ), (C.θ-1)/C.θ, 1)
 τ⁻¹(::Type{GumbelCopula},τ) =ifelse(τ == 1, Inf, 1/(1-τ))
-
-radial_dist(C::GumbelCopula) = AlphaStable(α = 1/C.θ, β = 1,scale = cos(π/(2C.θ))^C.θ, location = (C.θ == 1 ? 1 : 0))
+williamson_dist(C::GumbelCopula{d,T}) where {d,T} = WilliamsonFromFrailty(AlphaStable(α = 1/C.θ, β = 1,scale = cos(π/(2C.θ))^C.θ, location = (C.θ == 1 ? 1 : 0)), d)
 
 
 # S(α, β, γ , δ) denotes a stable distribution in

src/ArchimedeanCopulas/JoeCopula.jl

@@ -36,6 +36,6 @@ end
 ϕ(  C::JoeCopula,          t) = 1-(1-exp(-t))^(1/C.θ)
 ϕ⁻¹(C::JoeCopula,          t) = -log(1-(1-t)^C.θ)
 τ(C::JoeCopula) = 1 - 4sum(1/(k*(2+k*C.θ)*(C.θ*(k-1)+2)) for k in 1:1000) # 446 in R copula. 
-radial_dist(C::JoeCopula) = Sibuya(1/C.θ)
+williamson_dist(C::JoeCopula{d,T}) where {d,T} = WilliamsonFromFrailty(Sibuya(1/C.θ), d)
 
 

src/ArchimedeanCopulas/WilliamsonCopula.jl

@@ -0,0 +1,17 @@
+struct WilliamsonCopula{d,Tϕ,TX} <: ArchimedeanCopula{d}
+    ϕ::Tϕ
+    X::TX
+end
+function WilliamsonCopula(X::Distributions.UnivariateDistribution, d)
+    ϕ = WilliamsonTransforms.𝒲(X,d)
+    return WilliamsonCopula{d,typeof(ϕ),typeof(X)}(ϕ,X)
+end
+function WilliamsonCopula(ϕ::Function, d)
+    X = WilliamsonTransforms.𝒲₋₁(ϕ,d)
+    return WilliamsonCopula{d,typeof(ϕ),typeof(X)}(ϕ,X)
+end
+function WilliamsonCopula(ϕ::Function, X::Distributions.UnivariateDistribution, d)
+    return WilliamsonCopula{d,typeof(ϕ),typeof(X)}(ϕ,X)
+end
+williamson_dist(C::WilliamsonCopula) = C.X
+ϕ(C::WilliamsonCopula, t) = C.ϕ(t)

src/Copulas.jl

@@ -18,6 +18,7 @@ end Copulas
     import ForwardDiff
     import Cubature
     import MvNormalCDF
+    import WilliamsonTransforms
 
     # Standard copulas and stuff. 
     include("utils.jl")
@@ -50,6 +51,8 @@ end Copulas
     include("univariate_distributions/Sibuya.jl")
     include("univariate_distributions/Logarithmic.jl")
     include("univariate_distributions/AlphaStable.jl")
+    include("univariate_distributions/ClaytonWilliamsonDistribution.jl")
+    include("univariate_distributions/WilliamsonFromFrailty.jl")
 
     # Archimedean copulas
     include("ArchimedeanCopula.jl")
@@ -59,12 +62,13 @@ end Copulas
     include("ArchimedeanCopulas/GumbelCopula.jl")
     include("ArchimedeanCopulas/FrankCopula.jl")
     include("ArchimedeanCopulas/AMHCopula.jl")
+    include("ArchimedeanCopulas/WilliamsonCopula.jl")
     export IndependentCopula, 
            ClaytonCopula,
            JoeCopula,
            GumbelCopula,
            FrankCopula,
-           AMHCopula
-
+           AMHCopula,
+           WilliamsonCopula
 
 end

src/univariate_distributions/ClaytonWilliamsonDistribution.jl

@@ -0,0 +1,39 @@
+struct ClaytonWilliamsonDistribution{T<:Real,TI} <: Distributions.DiscreteUnivariateDistribution
+    θ::T
+    d::TI
+end
+function Distributions.cdf(D::ClaytonWilliamsonDistribution, x::Real)
+    θ = D.θ
+    d = D.d
+    if x < 0
+        return zero(x)
+    end
+    α = -1/θ
+    if θ < 0
+        if x >= α
+            return one(x)
+        end
+        rez = zero(x)
+        x_α = x/α
+        for k in 0:(d-1)
+            rez += SpecialFunctions.gamma(α+1)/SpecialFunctions.gamma(α-k+1)/SpecialFunctions.gamma(k+1) * (x_α)^k * (1 - x_α)^(α-k)
+        end
+        return 1-rez
+    elseif θ == 0
+        return exp(-x)
+    else
+        rez = zero(x)
+        for k in 0:(d-1)
+            pr = one(θ)
+            for j in 0:(k-1)
+                pr *= (1+j*θ)
+            end
+            rez += pr / SpecialFunctions.gamma(k+1) * x^k * (1 + θ * x)^(-(1/θ+k))
+        end
+        return 1-rez
+    end
+end
+function Distributions.rand(rng::Distributions.AbstractRNG, d::ClaytonWilliamsonDistribution)
+    u = rand(rng)
+    Roots.find_zero(x -> (Distributions.cdf(d,x) - u), (0.0, Inf))
+end

src/univariate_distributions/WilliamsonFromFrailty.jl

@@ -0,0 +1,11 @@
+struct WilliamsonFromFrailty{TF,d} <: Distributions.ContinuousUnivariateDistribution
+    frailty_dist::TF
+    function WilliamsonFromFrailty(frailty_dist,d)
+        return new{typeof(frailty_dist),d}(frailty_dist)
+    end
+end
+function Distributions.rand(rng::Distributions.AbstractRNG, D::WilliamsonFromFrailty{TF,d}) where {TF,d}
+    f = rand(rng,D.frailty_dist)
+    sy = sum(.-log.(rand(rng,d)))
+    return sy/f
+end