From 32e3e1c4bd0ff9ba0936352edd3d935174b28b5b Mon Sep 17 00:00:00 2001
From: Oskar Laverny <oskar.laverny@gmail.com>
Date: Thu, 19 Sep 2024 10:59:48 +0200
Subject: [PATCH] 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              | 17 +++++++++++++++++
 src/Copulas.jl                                  |  1 +
 .../UnivariateGenerator/GumbelGenerator.jl      |  9 +++++++--
 test/test_extreme_gumbel_density.jl             | 12 ++++++++++++
 4 files changed, 37 insertions(+), 2 deletions(-)
 create mode 100644 src/BivariateArchimedeanMethods.jl
 create mode 100644 test/test_extreme_gumbel_density.jl

diff --git a/src/BivariateArchimedeanMethods.jl b/src/BivariateArchimedeanMethods.jl
new file mode 100644
index 00000000..c13c9385
--- /dev/null
+++ b/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
diff --git a/src/Copulas.jl b/src/Copulas.jl
index 8f3fd192..9d887cf3 100644
--- a/src/Copulas.jl
+++ b/src/Copulas.jl
@@ -70,6 +70,7 @@ module Copulas
     
     # Archimedean copulas
     include("ArchimedeanCopula.jl")
+    include("BivariateArchimedeanMethods.jl")
     export ArchimedeanCopula,
            IndependentCopula,
            MCopula,
diff --git a/src/Generator/UnivariateGenerator/GumbelGenerator.jl b/src/Generator/UnivariateGenerator/GumbelGenerator.jl
index b797d751..2607f30c 100644
--- a/src/Generator/UnivariateGenerator/GumbelGenerator.jl
+++ b/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)
\ No newline at end of file
+williamson_dist(G::GumbelGenerator, d) = WilliamsonFromFrailty(AlphaStable(α = 1/G.θ, β = 1,scale = cos(π/(2G.θ))^G.θ, location = (G.θ == 1 ? 1 : 0)), d)
diff --git a/test/test_extreme_gumbel_density.jl b/test/test_extreme_gumbel_density.jl
new file mode 100644
index 00000000..fa3378d3
--- /dev/null
+++ b/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
\ No newline at end of file