Moving generators to their own structs

src/ArchimedeanCopula.jl

@@ -39,36 +39,36 @@ The Archimedean API is modular:
 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
+struct ArchimedeanCopula{d,TG} <: Copula{d}
+    G::TG
+    function ArchimedeanCopula(d::Int,G::Generator)
+        @assert d <= max_monotony(G) "The generator you provided is not d-monotonous according to its max_monotonicity property, and thus this copula does not exists."
+        return ArchimedeanCopula{d,typeof(G)}(G)
+    end
+    # three special cases: 
+    ArchimedeanCopula(d,::WGenerator) = WCopula(d)
+    ArchimedeanCopula(d,::IndependentGenerator) = IndependentCopula(d)
+    ArchimedeanCopula(d,::MGenerator) = MCopula(d)
+end
 function _cdf(C::CT,u) where {CT<:ArchimedeanCopula} 
     sum_ϕ⁻¹u = 0.0
     for us in u
-        sum_ϕ⁻¹u += ϕ⁻¹(C,us)
+        sum_ϕ⁻¹u += ϕ⁻¹(C.G,us)
     end
-    return ϕ(C,sum_ϕ⁻¹u)
-end
-ϕ⁽¹⁾(C::CT, t) where {CT<:ArchimedeanCopula} = ForwardDiff.derivative(x -> ϕ(C,x), t)
-function ϕ⁽ᵈ⁾(C::ArchimedeanCopula{d},t) where d
-    X = TaylorSeries.Taylor1(eltype(t),d)
-    taylor_expansion = ϕ(C,t+X)
-    coef = TaylorSeries.getcoeff(taylor_expansion,d) # gets the dth coef. 
-    der = coef * factorial(d) # gets the dth derivative of $\phi$ taken in t. 
-    return der
+    return ϕ(C.G,sum_ϕ⁻¹u)
 end
-function Distributions._logpdf(C::CT, u) where {CT<:ArchimedeanCopula}
-    d = length(C)
-    # @assert d == length(u) "Dimension mismatch"
+function Distributions._logpdf(C::ArchimedeanCopula{d,TG}, u) where {d,TG}
     if !all(0 .<= u .<= 1)
         return eltype(u)(-Inf)
     end
     sum_ϕ⁻¹u = 0.0
     sum_logϕ⁽¹⁾ϕ⁻¹u = 0.0
     for us in u
-        ϕ⁻¹u = ϕ⁻¹(C,us)
+        ϕ⁻¹u = ϕ⁻¹(C.G,us)
         sum_ϕ⁻¹u += ϕ⁻¹u
-        sum_logϕ⁽¹⁾ϕ⁻¹u += log(-ϕ⁽¹⁾(C,ϕ⁻¹u)) # log of negative here because ϕ⁽¹⁾ is necessarily negative
+        sum_logϕ⁽¹⁾ϕ⁻¹u += log(-ϕ⁽¹⁾(C.G,ϕ⁻¹u)) # log of negative here because ϕ⁽¹⁾ is necessarily negative
     end
-    numer = ϕ⁽ᵈ⁾(C, sum_ϕ⁻¹u)
+    numer = ϕ⁽ᵏ⁾(C.G, d, sum_ϕ⁻¹u)
     dimension_sign = iseven(d) ? 1.0 : -1.0 #need this for log since (-1.0)ᵈ ϕ⁽ᵈ⁾ ≥ 0.0
 
 
@@ -83,35 +83,31 @@ function Distributions._logpdf(C::CT, u) where {CT<:ArchimedeanCopula}
         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 = Roots.find_zero(t -> ϕ(C,t) - x, (0.0, Inf))
-
-williamson_dist(C::ArchimedeanCopula{d}) where d = WilliamsonTransforms.𝒲₋₁(t -> ϕ(C,t),d)
+_W(C::ArchimedeanCopula{d,TG}) where {d,TG} = williamson_dist(C.G,d)
 function τ(C::ArchimedeanCopula)  
     @show C
-    return 4*Distributions.expectation(r -> ϕ(C,r), williamson_dist(C)) - 1
+    return 4*Distributions.expectation(r -> ϕ(C,r), _W(C)) - 1
 end
 
 function _archi_rand!(rng,C::ArchimedeanCopula{d},R,x) where d
     # x is assumed to already be random exponentials produced by Random.randexp
     r = rand(rng,R)
     sx = sum(x)
     for i in 1:d
-        x[i] = ϕ(C,r * x[i]/sx)
+        x[i] = ϕ(C.G,r * x[i]/sx)
     end
 end
 
 function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, x::AbstractVector{T}) where {T<:Real, CT<:ArchimedeanCopula}
     # By default, we use the williamson sampling. 
     Random.randexp!(rng,x)
-    _archi_rand!(rng,C,williamson_dist(C),x)
+    _archi_rand!(rng,C.G,_W(C),x)
     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: 
     Random.randexp!(rng,A)
-    R = williamson_dist(C)
+    R = _W(C)
     for i in 1:size(A,2)
         _archi_rand!(rng,C,R,view(A,:,i))
     end

src/ArchimedeanCopulas/AMHCopula.jl

@@ -17,27 +17,11 @@ The [AMH](https://en.wikipedia.org/wiki/Copula_(probability_theory)#Most_importa
 It has a few special cases: 
 - When θ = 0, it is the IndependentCopula
 """
-struct AMHCopula{d,T} <: ArchimedeanCopula{d}
-    θ::T
-    function AMHCopula(d,θ)
-        if (θ < -1) || (θ >= 1)
-            throw(ArgumentError("Theta must be in [-1,1)"))
-        elseif θ == 0
-            return IndependentCopula(d)
-        else
-            return new{d,typeof(θ)}(θ)
-        end
-    end
-end
-ϕ(  C::AMHCopula,t) = (1-C.θ)/(exp(t)-C.θ)
-ϕ⁻¹(  C::AMHCopula,t) = log(C.θ + (1-C.θ)/t)
-
-τ(C::AMHCopula) = _amh_tau_f(C.θ) # no closed form inverse...
+const AMHCopula{d,T} = ArchimedeanCopula{d,AMHGenerator{T}}
+AMHCopula(d,θ) = ArchimedeanCopula(d,AMHGenerator(θ))
 
 _amh_tau_f(θ) = θ == 1 ? 1/3 : 1 - 2(θ+(1-θ)^2*log1p(-θ))/(3θ^2)
-
-# if θ = -1, we obtain (5 -8*log(2))/3 
-
+τ(C::AMHCopula) = _amh_tau_f(C.θ) # no closed form inverse...
 function τ⁻¹(::Type{AMHCopula},τ)
     if τ == zero(τ)
         return τ
@@ -51,7 +35,4 @@ function τ⁻¹(::Type{AMHCopula},τ)
         return -one(τ)
     end
     return Roots.find_zero(θ -> _amh_tau_f(θ) - τ, (-one(τ), one(τ)))
-end
-williamson_dist(C::AMHCopula{d,T}) where {d,T} = C.θ >= 0 ? WilliamsonFromFrailty(1 + Distributions.Geometric(1-C.θ),d) : WilliamsonTransforms.𝒲₋₁(t -> ϕ(C,t),d)
-
-
+end

src/ArchimedeanCopulas/ClaytonCopula.jl

@@ -19,24 +19,7 @@ It has a few special cases:
 - When θ = 0, it is the IndependentCopula
 - When θ = ∞, is is the MCopula (Upper Frechet-Hoeffding bound)
 """
-struct ClaytonCopula{d,T} <: ArchimedeanCopula{d}
-    θ::T
-    function ClaytonCopula(d,θ)
-        if θ < -1/(d-1)
-            throw(ArgumentError("Theta must be greater than -1/(d-1)"))
-        elseif θ == -1/(d-1)
-            return WCopula(d)
-        elseif θ == 0
-            return IndependentCopula(d)
-        elseif θ == Inf
-            return MCopula(d)
-        else
-            return new{d,typeof(θ)}(θ)
-        end
-    end
-end
-ϕ(  C::ClaytonCopula,      t) = max(1+C.θ*t,zero(t))^(-1/C.θ)
-ϕ⁻¹(C::ClaytonCopula,      t) = (t^(-C.θ)-1)/C.θ
+const ClaytonCopula{d,T} = ArchimedeanCopula{d,ClaytonGenerator{T}}
+ClaytonCopula(d,θ) = ArchimedeanCopula(d,ClaytonGenerator(θ))
 τ(C::ClaytonCopula) = ifelse(isfinite(C.θ), C.θ/(C.θ+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.θ,C.θ),d) : ClaytonWilliamsonDistribution(C.θ,d)
+τ⁻¹(::Type{ClaytonCopula},τ) = ifelse(τ == 1,Inf,2τ/(1-τ))

src/ArchimedeanCopulas/FrankCopula.jl

@@ -19,40 +19,13 @@ It has a few special cases:
 - When θ = 1, it is the IndependentCopula
 - When θ = ∞, is is the MCopula (Upper Frechet-Hoeffding bound)
 """
-struct FrankCopula{d,T} <: ArchimedeanCopula{d}
-    θ::T
-    function FrankCopula(d,θ)
-        if d > 2 && θ < 0
-            throw(ArgumentError("Negatively dependent Frank copulas cannot exists in dimensions > 2"))
-        end
-        if θ == -Inf
-            return WCopula(d)
-        elseif θ == 0
-            return IndependentCopula(d)
-        elseif θ == Inf
-            return MCopula(d)
-        else
-            return new{d,typeof(θ)}(θ)
-        end
-    end
-end
-ϕ(  C::FrankCopula,       t) = C.θ > 0 ? -LogExpFunctions.log1mexp(LogExpFunctions.log1mexp(-C.θ)-t)/C.θ : -log1p(exp(-t) * expm1(-C.θ))/C.θ
-ϕ⁻¹(C::FrankCopula,       t) = C.θ > 0 ? LogExpFunctions.log1mexp(-C.θ) - LogExpFunctions.log1mexp(-t*C.θ) : -log(expm1(-t*C.θ)/expm1(-C.θ))
-
-# A bit of type piracy but should be OK : 
-# LogExpFunctions.log1mexp(t::TaylorSeries.Taylor1) = log(-expm1(t))
-# Avoid type piracy by defiing it myself: 
-ϕ(  C::FrankCopula,       t::TaylorSeries.Taylor1) = C.θ > 0 ? -log(-expm1(LogExpFunctions.log1mexp(-C.θ)-t))/C.θ : -log1p(exp(-t) * expm1(-C.θ))/C.θ
+const FrankCopula{d,T} = ArchimedeanCopula{d,FrankGenerator{T}}
+FrankCopula(d,θ) = ArchimedeanCopula(d,FrankGenerator(θ))
 
 D₁ = GSL.sf_debye_1 # sadly, this is C code.
 # could be replaced by : 
 # using QuadGK
 # D₁(x) = quadgk(t -> t/(exp(t)-1), 0, x)[1]/x
-# to make it more general. but once gain, it requires changing the integrator at each evlauation, 
-# which is problematic. 
-# Better option is to try to include this function into SpecialFunctions.jl. 
-
-
 τ(C::FrankCopula) = 1+4(D₁(C.θ)-1)/C.θ
 function τ⁻¹(::Type{FrankCopula},τ)
     if τ == zero(τ)
@@ -63,7 +36,4 @@ function τ⁻¹(::Type{FrankCopula},τ)
     end
     x₀ = (1-τ)/4
     return Roots.fzero(x -> (1-D₁(x))/x - x₀, 1e-4, Inf)
-end
-
-williamson_dist(C::FrankCopula{d,T}) where {d,T} = C.θ > 0 ?  WilliamsonFromFrailty(Logarithmic(-C.θ), d) : WilliamsonTransforms.𝒲₋₁(t -> ϕ(C,t),d)
-
+end

src/ArchimedeanCopulas/IndependentCopula.jl

@@ -34,4 +34,11 @@ end
 
 function Distributions._rand!(rng::Distributions.AbstractRNG, C::IndependentCopula{d}, x::AbstractVector{T}) where {T<:Real,d}
     Random.rand!(rng,x)
-end
+end
+
+
+######## maybe this should NOT be an archimedean copula, but rather stay with the W / Pi / M copulas in the miscelaneous stuff ? 
+######## Or on the other way around,maybe W Pi AND M should be in the archimdean class ? even if M is not really archimedean... 
+######## For M, the genreatormight have a weird head.. 
+
+

src/Generator.jl

@@ -0,0 +1,32 @@
+abstract type Generator end
+
+max_monotony(G::Generator) = throw("This generator does not have a defined max monotony")
+ϕ(   G::Generator, t) = throw("This generator has not been defined correctly.")
+ϕ⁻¹( G::Generator, x) = Roots.find_zero(t -> ϕ(G,t) - x, (0.0, Inf))
+ϕ⁽¹⁾(G::Generator, t) = ForwardDiff.derivative(x -> ϕ(G,x), t)
+function ϕ⁽ᵏ⁾(G::Generator, k, t)
+    X = TaylorSeries.Taylor1(eltype(t),k)
+    taylor_expansion = ϕ(C,t+X)
+    coef = TaylorSeries.getcoeff(taylor_expansion,k) 
+    der = coef * factorial(d)
+    return der
+end
+williamson_dist(G::Generator, d) = WilliamsonTransforms.𝒲₋₁(t -> ϕ(G,t),d)
+
+
+# Maybe the three following generators do not need to be that much special case ? 
+# Maybe they do. 
+
+
+struct WGenerator <: Generator end
+# max_monotony(G::Wgenerator) = 2
+# ϕ(G::WGenerator, t) = 
+
+struct IndependentGenerator <: Generator end
+# max_monotony(G::IndependentGenerator) = Inf
+# ϕ(   G::IndependentGenerator, t) = exp(-t)
+
+struct MGenerator <: Generator end
+# max_monotony(G::MGenerator) = Inf
+# ϕ(   G::MGenerator, t) = throw(ArgumentError("The MGenerator should never be called."))
+

src/Generators/AMHGenerator.jl

@@ -0,0 +1,18 @@
+struct AMHGenerator{T} <: Generator
+    θ::T
+    function AMHGenerator(θ)
+        if (θ < -1) || (θ >= 1)
+            throw(ArgumentError("Theta must be in [-1,1)"))
+        elseif θ == 0
+            return IndependentGenrator()
+        else
+            return new{d,typeof(θ)}(θ)
+        end
+    end
+end
+max_monotony(G::AMHGenerator) = Inf
+ϕ(  G::AMHGenerator, t) = (1-G.θ)/(exp(t)-G.θ)
+ϕ⁻¹(G::AMHGenerator, t) = log(G.θ + (1-G.θ)/t)
+# ϕ⁽¹⁾(G::AMHGenerator, t) =  First derivative of ϕ
+# ϕ⁽ᵏ⁾(G::AMHGenerator, k, t) = kth derivative of ϕ
+williamson_dist(G::AMHGenerator, d) = G.θ >= 0 ? WilliamsonFromFrailty(1 + Distributions.Geometric(1-G.θ),d) : WilliamsonTransforms.𝒲₋₁(t -> ϕ(G,t),d)

src/Generators/ClaytonGenerator.jl

@@ -0,0 +1,22 @@
+struct ClaytonGenerator{T} <: Generator
+    θ::T
+    function ClaytonGenerator(θ)
+        if θ < -1
+            throw(ArgumentError("Theta must be greater than -1"))
+        elseif θ == -1
+            return WGenerator()
+        elseif θ == 0
+            return IndependentGenrator()
+        elseif θ == Inf
+            return MGenerator()
+        else
+            return new{typeof(θ)}(θ)
+        end
+    end
+end
+max_monotony(G::ClaytonGenerator) = Int(floor(1 - 1/G.θ))
+ϕ(  G::ClaytonGenerator, t) = max(1+G.θ*t,zero(t))^(-1/G.θ)
+ϕ⁻¹(G::ClaytonGenerator, t) = (t^(-G.θ)-1)/G.θ
+# ϕ⁽¹⁾(G::ClaytonGenerator, t) =  First derivative of ϕ
+# ϕ⁽ᵏ⁾(G::ClaytonGenerator, k, t) = kth derivative of ϕ
+williamson_dist(G::ClaytonGenerator, d) = G.θ >= 0 ? WilliamsonFromFrailty(Distributions.Gamma(1/C.θ,C.θ),d) : ClaytonWilliamsonDistribution(C.θ,d)

src/Generators/FrankGenerator.jl

@@ -0,0 +1,21 @@
+struct FrankGenerator{T} <: Generator
+    θ::T
+    function FrankGenerator(θ)
+        if θ == -Inf
+            return WGenerator()
+        elseif θ == 0
+            return IndependentGenrator()
+        elseif θ == Inf
+            return MGenerator()
+        else
+            return new{d,typeof(θ)}(θ)
+        end
+    end
+end
+max_monotony(G::FrankGenerator) = G.θ < 0 ? 2 : Inf
+ϕ(  G::FrankGenerator, t) = G.θ > 0 ? -LogExpFunctions.log1mexp(LogExpFunctions.log1mexp(-G.θ)-t)/G.θ : -log1p(exp(-t) * expm1(-G.θ))/G.θ
+ϕ(  G::FrankGenerator, t::TaylorSeries.Taylor1) = G.θ > 0 ? -log(-expm1(LogExpFunctions.log1mexp(-G.θ)-t))/G.θ : -log1p(exp(-t) * expm1(-G.θ))/G.θ
+ϕ⁻¹(G::FrankGenerator, t) = G.θ > 0 ? LogExpFunctions.log1mexp(-G.θ) - LogExpFunctions.log1mexp(-t*G.θ) : -log(expm1(-t*G.θ)/expm1(-G.θ))
+# ϕ⁽¹⁾(G::FrankGenerator, t) =  First derivative of ϕ
+# ϕ⁽ᵏ⁾(G::FrankGenerator, k, t) = kth derivative of ϕ
+williamson_dist(G::FrankGenerator, d) = G.θ > 0 ?  WilliamsonFromFrailty(Logarithmic(-G.θ), d) : WilliamsonTransforms.𝒲₋₁(t -> ϕ(G,t),d)