Fix galambox (#220)

* modified:   src/ExtremeValueCopulas/GalambosCopula.jl
	modified:   src/UnivariateDistribution/ExtremeDist.jl
	modified:   test/Extreme_value_test.jl

* fix constructor

* removed doubled test

* add test

* refactor extreme values tests

* typo

---------

Co-authored-by: Santymax98 <[email protected]>

src/ExtremeValueCopulas/GalambosCopula.jl

@@ -26,6 +26,9 @@ References:
 struct GalambosCopula{P} <: ExtremeValueCopula{P}
     θ::P  # Copula parameter
     function GalambosCopula(θ)
+        if θ > 19.5
+            @warn "The parameter θ = $(θ) is large, which may lead to numerical instability in any functions. Consider regularizing the input."
+        end
         if θ < 0
             throw(ArgumentError("Theta must be >= 0"))
         elseif θ == 0
@@ -38,4 +41,8 @@ struct GalambosCopula{P} <: ExtremeValueCopula{P}
     end
 end
 
-A(C::GalambosCopula, t::Real) = 1 - (t^(-C.θ) + (1 - t)^(-C.θ))^(-1/C.θ)
+A(C::GalambosCopula, t::Real) = 1 - (t^(-C.θ) + (1 - t)^(-C.θ))^(-1/C.θ)
+# This auxiliary function helps determine if we need binary search or not in the generation of random samples
+function needs_binary_search(C::GalambosCopula)
+    return C.θ > 19.5
+end

src/UnivariateDistribution/ExtremeDist.jl

@@ -6,14 +6,11 @@ struct ExtremeDist{C} <: Distributions.ContinuousUnivariateDistribution
 end
 
 function Distributions.cdf(d::ExtremeDist, z)
-    if z < 0
+    if z < 0 || z > 1
         return 0.0
-    elseif z > 1
-        return 1.0
-    else
-        copula = d.G
-        return z + z * (1 - z) * (dA(copula, z) / A(copula, z))
     end
+    copula = d.G
+    return z + z * (1 - z) * (dA(copula, z)/A(copula, z))
 end
 
 function _pdf(d::ExtremeDist, z)
@@ -33,7 +30,34 @@ function Distributions.quantile(d::ExtremeDist, p)
         throw(ArgumentError("p must be between 0 and 1"))
     end
     cdf_func(x) = Distributions.cdf(d, x) - p
-    return Roots.find_zero(cdf_func, (eps(), 1-eps()), Roots.Brent())
+    copula = d.G
+
+    # Automatically decide whether to use binary search or Brent
+    if hasmethod(needs_binary_search, (typeof(copula),)) && needs_binary_search(copula)
+        # Use binary search for copulas with large parameters     
+        lower_bound = eps()
+        upper_bound = 1.0 - eps()
+        mid_point = (lower_bound + upper_bound) / 2
+
+        while upper_bound - lower_bound > 1e-6  # Accuracy threshold
+            mid_value = cdf_func(mid_point)
+
+            if abs(mid_value) < 1e-6
+                return mid_point
+            elseif mid_value > 0
+                upper_bound = mid_point
+            else
+                lower_bound = mid_point
+            end
+
+            mid_point = (lower_bound + upper_bound) / 2
+        end
+
+        return mid_point
+    else
+        # Use Brent for other copulations or if there are no problems
+        return Roots.find_zero(cdf_func, (eps(), 1.0 - eps()), Roots.Brent())
+    end
 end
 
 # Generate random samples from the radial distribution using the quantile function