Frank fix + coherence tests on Archimedeans + marginals tests on all copulas

Project.toml

@@ -9,6 +9,7 @@ Cubature = "667455a9-e2ce-5579-9412-b964f529a492"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 GSL = "92c85e6c-cbff-5e0c-80f7-495c94daaecd"
+LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 MvNormalCDF = "37188c8d-bc69-4638-b057-733e744175ec"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
@@ -18,17 +19,24 @@ TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea"
 WilliamsonTransforms = "48feb556-9bdd-43a2-8e10-96100ec25e22"
 
 [compat]
+Aqua = "v0.8"
 Combinatorics = "1"
 Cubature = "1.5"
 Distributions = "0.25"
 ForwardDiff = "0.10"
 GSL = "1"
+HypothesisTests = "v0.11"
+InteractiveUtils = "1.6"
+LinearAlgebra = "1.6"
+LogExpFunctions = "0.3"
 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"
+Test = "1.6"
+TestItemRunner = "v0.2"
 WilliamsonTransforms = "0.1"
 julia = "1.6"
 

src/ArchimedeanCopulas/ClaytonCopula.jl

@@ -39,4 +39,4 @@ end
 ϕ⁻¹(C::ClaytonCopula,      t) = (t^(-C.θ)-1)/C.θ
 τ(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.θ,1),d) : ClaytonWilliamsonDistribution(C.θ,d)
+williamson_dist(C::ClaytonCopula{d,T}) where {d,T} = C.θ >= 0 ? WilliamsonFromFrailty(Distributions.Gamma(1/C.θ,C.θ),d) : ClaytonWilliamsonDistribution(C.θ,d)

src/ArchimedeanCopulas/FrankCopula.jl

@@ -22,6 +22,9 @@
 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
@@ -33,8 +36,13 @@
         end
     end
 end
-ϕ(  C::FrankCopula,       t) = -log(1+exp(-t)*(exp(-C.θ)-1))/C.θ
-ϕ⁻¹(C::FrankCopula,       t) = -log((exp(-t*C.θ)-1)/(exp(-C.θ)-1))
+ϕ(  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.θ
 
 D₁ = GSL.sf_debye_1 # sadly, this is C code.
 # could be replaced by : 
@@ -54,6 +62,5 @@
     return Roots.fzero(x -> (1-D₁(x))/x - x₀, 1e-4, Inf)
 end
 
-williamson_dist(C::FrankCopula{d,T}) where {d,T} = WilliamsonFromFrailty(Logarithmic(1-exp(-C.θ)), d)
-
+williamson_dist(C::FrankCopula{d,T}) where {d,T} = C.θ > 0 ?  WilliamsonFromFrailty(Logarithmic(-C.θ), d) : WilliamsonTransforms.𝒲₋₁(t -> ϕ(C,t),d)
 

src/ArchimedeanCopulas/JoeCopula.jl

@@ -32,8 +32,8 @@ struct JoeCopula{d,T} <: ArchimedeanCopula{d}
         end
     end
 end
-ϕ(  C::JoeCopula,          t) = 1-(1-exp(-t))^(1/C.θ)
-ϕ⁻¹(C::JoeCopula,          t) = -log(1-(1-t)^C.θ)
+ϕ(  C::JoeCopula,          t) = 1-(-expm1(-t))^(1/C.θ)
+ϕ⁻¹(C::JoeCopula,          t) = -log1p(-(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. 
 williamson_dist(C::JoeCopula{d,T}) where {d,T} = WilliamsonFromFrailty(Sibuya(1/C.θ), d)
 

src/Copulas.jl

@@ -13,6 +13,7 @@ module Copulas
     import MvNormalCDF
     import WilliamsonTransforms
     import Combinatorics
+    import LogExpFunctions
 
     # Standard copulas and stuff. 
     include("utils.jl")

src/univariate_distributions/AlphaStable.jl

@@ -29,7 +29,7 @@ function Distributions.rand(rng::Distributions.AbstractRNG, d::AlphaStable{T}) w
 # McCulloch's MATLAB implementation (1996) served as a reference in developing this code.
 # """
     α=d.α; β=d.β; sc=d.scale; loc=d.location
-    (α < 0.1 || α > 2) && throw(DomainError(α, "α must be in the range 0.1 to 2"))
+    (α <= 0 || α > 2) && throw(DomainError(α, "α must be in the range 0.1 to 2"))
     abs(β) > 1 && throw(DomainError(β, "β must be in the range -1 to 1"))
     # added eps(T) to prevent DomainError: x ^ y where x < 0
     ϕ = (rand(rng, T) - T(0.5)) * π * (one(T) - eps(T))

src/univariate_distributions/Logarithmic.jl

@@ -1,25 +1,54 @@
 # Corresponds to https://en.wikipedia.org/wiki/Logarithmic_distribution
 struct Logarithmic{T<:Real} <: Distributions.DiscreteUnivariateDistribution
-    p::T
-    function Logarithmic(p::T) where {T <: Real}
-        new{T}(p)
+    α::T # in (0,1), the weight of the logarithmic distribution.
+    h::T # = ln(1-α), so in (-Inf,0)
+    function Logarithmic(h::T) where {T <: Real}
+        Tf = promote_type(T,Float64)
+        α = -expm1(h)
+        return new{Tf}(Tf(α), Tf(h))
     end
 end
-function Distributions.logpdf(d::Logarithmic, x::Real)
-    insupport(d, x) ? x*log(d.p) - log(x) - log(-log(1-d.p)) : log(zero(d.p))
+Base.eltype(::Logarithmic{T}) where T = promote_type(T,Float64)
+function Distributions.logpdf(d::Logarithmic{T}, x::Real) where T
+    insupport(d, x) ? x*log1p(-d.α) - log(x) - log(-log(d.α)) : log(zero(T))
 end
-function Distributions.rand(rng::Distributions.AbstractRNG, d::Logarithmic)
-    # Sample a Log(p) distribution with the algorithm "LK" of Kemp (1981).
-    u = rand(rng)
-    if u > d.p 
-        return 1
+function Distributions.rand(rng::Distributions.AbstractRNG, d::Logarithmic{T}) where T
+    # Sample a Log(p) distribution with the algorithms "LK" and "LS" of Kemp (1981).
+    if d.h > -3
+        # Version "LS"
+        t = -d.α / log1p(-d.α)
+        u = rand(rng)
+        p = t
+        x = 1
+        while u > p
+            u = u - p
+            x = x + 1
+            p = p * d.α * (x-1) / x
+        end
+        return x
+    else
+        # Version "LK"
+        u = rand(rng)
+        if u > d.α
+            return one(T)
+        else
+            v = rand(rng)
+            h2 = d.h * v
+            r = -expm1(h2)
+            if u < r*r 
+                lr = LogExpFunctions.log1mexp(h2)
+                if iszero(lr)
+                    return T(Inf)
+                else
+                    return floor(1 + log(u)/lr)
+                end
+            else
+                if u > r 
+                    return one(T)
+                else
+                    return T(2)
+                end
+            end
+        end
     end
-    q = 1 - (1-d.p)^rand(rng)
-    if u < q*q
-        return floor(1+log(u)/log(q))
-    end
-    if u < q
-        return 1
-    end
-    return 2
-end
+end

test/archimedean_tests.jl

@@ -42,8 +42,9 @@ end
 end
 @testitem "FrankCopula - Test sampling all cases" begin
     rand(FrankCopula(2,-Inf),10)
+    rand(FrankCopula(2,log(rand())),10)
     for d in 2:10
-        for θ ∈ [log(rand()),0.0,rand(),1.0,-log(rand()), Inf]
+        for θ ∈ [0.0,rand(),1.0,-log(rand()), Inf]
             rand(FrankCopula(d,θ),10)
         end
     end

test/coherence_archimedeans.jl

@@ -1,17 +1,24 @@
 
 @testitem "Test of coherence for archimedeans" begin
-    using HypothesisTests, Distributions
+    using HypothesisTests, Distributions, Random
+    Random.seed!(123)
     cops = (
         # true represent the fact that cdf(williamson_dist(C),x) is defined or not. 
         (AMHCopula(3,0.6), true),
         (AMHCopula(4,-0.3), true),
         (ClaytonCopula(2,-0.7), true),
         (ClaytonCopula(3,-0.1), true),
         (ClaytonCopula(4,7), true),
-        (FrankCopula(3,-12), false),
+        (FrankCopula(2,-5), false),
+        (FrankCopula(3,12), false),
         (FrankCopula(4,6), false),
+        (FrankCopula(4,30), false),
+        (FrankCopula(4,37), false),
+        (FrankCopula(4,150), false),
         (JoeCopula(3,7), false),
         (GumbelCopula(4,7), false),
+        (GumbelCopula(4,20), false),
+        (GumbelCopula(4,100), false),
         (GumbelBarnettCopula(3,0.7),true),
         (InvGaussianCopula(4,0.05),true),
         (InvGaussianCopula(3,8),true)
@@ -20,16 +27,16 @@
     spl = rand(n)
     spl2 = rand(n)
     for (C,will_dist_has_a_cdf_implemented) in cops
-        spl .= dropdims(sum(Copulas.ϕ.(Ref(C),rand(C,n)),dims=1),dims=1)
+        spl .= dropdims(sum(Copulas.ϕ⁻¹.(Ref(C),rand(C,n)),dims=1),dims=1)
         will_dist = Copulas.williamson_dist(C)
         if will_dist_has_a_cdf_implemented
-            pval = pvalue(ExactOneSampleKSTest(spl, will_dist))
-            @test pval < 0.05
+            pval = pvalue(ExactOneSampleKSTest(spl, will_dist),tail=:right)
+            @test pval > 0.01
         end
         # even without a cdf we can still test approximately:
         spl2 .= rand(will_dist,n)
-        pval2 = pvalue(ApproximateTwoSampleKSTest(spl,spl2))
-        @test pval2 < 0.05
+        pval2 = pvalue(ApproximateTwoSampleKSTest(spl,spl2),tail=:right)
+        @test pval2 > 0.01
     end
 
 end

test/margins_uniformity.jl

@@ -0,0 +1,41 @@
+@testitem "Test samples have uniform maginals in [0,1]" begin
+    using HypothesisTests, Distributions, Random
+    Random.seed!(1234)
+    cops = (
+        # true represent the fact that cdf(williamson_dist(C),x) is defined or not. 
+        AMHCopula(3,0.6),
+        AMHCopula(4,-0.3),
+        ClaytonCopula(2,-0.7),
+        ClaytonCopula(3,-0.1),
+        ClaytonCopula(4,7),
+        FrankCopula(2,-5),
+        FrankCopula(3,12),
+        FrankCopula(4,6),
+        FrankCopula(4,150),
+        JoeCopula(3,7),
+        GumbelCopula(4,7),
+        GumbelCopula(4,20),
+        GumbelCopula(4,100),
+        GumbelBarnettCopula(3,0.7),
+        InvGaussianCopula(4,0.05),
+        InvGaussianCopula(3,8),
+        GaussianCopula([1 0.5; 0.5 1]),
+        TCopula(4, [1 0.5; 0.5 1]),
+        FGMCopula(2,1),
+        MCopula(4),
+        PlackettCopula(2.0),
+        # Others ? Yes probably others too ! 
+    )
+    n = 1000
+    U = Uniform(0,1)
+    for C in cops
+        nfail = 0
+        d = length(C)
+        @show C
+        spl = rand(C,n)
+        @assert all(0 <= x <= 1 for x in spl)
+        for i in 1:d
+            @test pvalue(ApproximateOneSampleKSTest(spl[i,:], U),tail=:right) > 0.01 # quite weak but enough at these samples sizes to detect really bad behaviors.
+        end
+    end
+end

test/test_extreme_frank_density.jl

@@ -0,0 +1,6 @@
+@testitem "Extreme frank density test" begin
+    using Distributions
+    F = FrankCopula(2, 60)
+    den = pdf(F,[0.70, 0.66])
+    @test true
+end