Remove overflows for bivariate gumbel density (#218)

* implement bivariate logdensity for gumbel to remove overflows.
Fixes #215

* add also the cdf

* add the cdf too

* move the code around

* nailed it !

* change comment

* Add a specific file for bivariate archimedean methods

src/BivariateArchimedeanMethods.jl

@@ -0,0 +1,17 @@
+# 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}
+    !all(0 .<= u .<= 1) && return eltype(u)(-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

src/Copulas.jl

@@ -70,6 +70,7 @@ module Copulas
 
     # Archimedean copulas
     include("ArchimedeanCopula.jl")
+    include("BivariateArchimedeanMethods.jl")
     export ArchimedeanCopula,
            IndependentCopula,
            MCopula,

src/Generator/UnivariateGenerator/GumbelGenerator.jl

@@ -39,7 +39,12 @@ end
 max_monotony(G::GumbelGenerator) = Inf
 ϕ(  G::GumbelGenerator, t) = exp(-t^(1/G.θ))
 ϕ⁻¹(G::GumbelGenerator, t) = (-log(t))^G.θ
-# ϕ⁽¹⁾(G::GumbelGenerator, t) =  First derivative of ϕ
+function ϕ⁽¹⁾(G::GumbelGenerator, t)
+    # first derivative of ϕ
+    a = 1/G.θ
+    tam1 = t^(a-1)
+    return - a * tam1 * exp(-tam1*t)
+end
 # ϕ⁽ᵏ⁾(G::GumbelGenerator, k, t) = kth derivative of ϕ
 τ(G::GumbelGenerator) = ifelse(isfinite(G.θ), (G.θ-1)/G.θ, 1)
 function τ⁻¹(::Type{T},τ) where T<:GumbelGenerator 
@@ -54,4 +59,4 @@ function τ⁻¹(::Type{T},τ) where T<:GumbelGenerator
         return θ
     end
 end
-williamson_dist(G::GumbelGenerator, d) = WilliamsonFromFrailty(AlphaStable(α = 1/G.θ, β = 1,scale = cos(π/(2G.θ))^G.θ, location = (G.θ == 1 ? 1 : 0)), d)
+williamson_dist(G::GumbelGenerator, d) = WilliamsonFromFrailty(AlphaStable(α = 1/G.θ, β = 1,scale = cos(π/(2G.θ))^G.θ, location = (G.θ == 1 ? 1 : 0)), d)

test/test_extreme_gumbel_density.jl

@@ -0,0 +1,12 @@
+@testitem "Extreme gumbel density test" begin
+    using Distributions
+    G = GumbelCopula(2, 244.3966206319112)
+    u = [0.969034207932297,0.9638122014293545]
+    den = pdf(G, u)
+    @test true
+
+    G = GumbelCopula(2, 1544.3966206319112)
+    u = [0.969034207932297,0.9638122014293545]
+    den = pdf(G, u)
+    @test true
+end