[Bug] RafteryCopula fixes. (#137)

* Fixing the pdf and clearing out the code.
Fixes #98

* Small fix for boundary conditions

* Add final test and fix previous values

* Add manually computed values

---------

Co-authored-by: Santymax98

src/Copula.jl

@@ -11,18 +11,18 @@ Base.length(::Copula{d}) where d = d
 # Base.eltype 
 function Distributions.cdf(C::Copula{d},u::VT) where {d,VT<:AbstractVector}
     length(u) != d && throw(ArgumentError("Dimension mismatch between copula and input vector"))
+    if any(iszero,u)
+        return zero(u[1])
+    elseif all(isone,u)
+        return one(u[1])
+    end
     return _cdf(C,u)
 end
 function Distributions.cdf(C::Copula{d},A::AbstractMatrix) where d
     size(A,1) != d && throw(ArgumentError("Dimension mismatch between copula and input vector"))
-    return [_cdf(C,u) for u in eachcol(A)]
+    return [Distributions.cdf(C,u) for u in eachcol(A)]
 end
 function _cdf(C::CT,u) where {CT<:Copula}
-    if any(iszero.(u))
-        return zero(u[1])
-    elseif all(isone.(u))
-        return one(u[1])
-    end
     f(x) = Distributions.pdf(C,x)
     z = zeros(eltype(u),length(C))
     return Cubature.hcubature(f,z,u,reltol=sqrt(eps()))[1]

src/MiscellaneousCopulas/RafteryCopula.jl

@@ -8,7 +8,7 @@ Constructor
 
     RafteryCopula(d, θ)
 
-The Multivariate Raftery Copula of dimension d is arameterized by ``\\theta \\in [0,1]`` 
+The Multivariate Raftery Copula of dimension d is parameterized by ``\\theta \\in [0,1]`` 
 
 ```math
 C_{\\theta}(\\mathbf{u}) = u_{(1)} + \\frac{(1 - \\theta)(1 - d)}{1 - \\theta - d} \\left(\\prod_{j=1}^{d} u_j\\right)^{\\frac{1}{1-\\theta}} - \\sum_{i=2}^{d} \\frac{\\theta(1-\\theta)}{(1-\\theta-i)(2-\\theta-i)} \\left(\\prod_{j=1}^{i-1}u_{(j)}\\right)^{\\frac{1}{1-\\theta}}u_{(i)}^{\\frac{2-\\theta-i}{1-\\theta}}
@@ -39,58 +39,37 @@ struct RafteryCopula{d, P} <: Copula{d}
     end
 end
 Base.eltype(R::RafteryCopula) = eltype(R.θ)
-
 function _cdf(R::RafteryCopula{d,P}, u) where {d,P}
-
-    if any(iszero,u)
-        return zero(u[1])
-    end
-    if all(isone,u)
-        return one(u[1])
-    end
-
     # Order the vector u
     u_ordered = sort(u)
-
     term1 = u_ordered[1]
     term2 = (1 - R.θ) * (1 - d) / (1 - R.θ - d) * prod(u).^(1/(1 - R.θ))
-
     term3 = 0.0
     for i in 2:d
         prod_prev = prod(u_ordered[1:i-1])
-        term3 += R.θ * (1 - R.θ) / ((1 - R.θ - i) * (2 - R.θ - i)) * prod_prev^(1/(1 - R.θ)) * u_ordered[i]^((2 - R.θ - i) / (1 - R.θ))
+        term3_part = R.θ * (1 - R.θ) / ((1 - R.θ - i) * (2 - R.θ - i)) * prod_prev^(1/(1 - R.θ)) * u_ordered[i]^((2 - R.θ - i) / (1 - R.θ))
+        term3 += term3_part
     end
-    # Combine the terms to get the cumulative distribution function
-    cdf_value = term1 + term2 - term3
-
-    return cdf_value
+    return term1 + term2 - term3
 end
 function Distributions._logpdf(R::RafteryCopula{d,P}, u::Vector{T}) where {d,P,T}
     # Order the vector u
     u_ordered = sort(u)
-
-    term_denominator = (1 - R.θ)^(d - 1) * (1 - R.θ - d)
-    term_numerator = 1 - d - R.θ * u_ordered[d]^((1 - R.θ - d) / (1 - R.θ))
-    term_product = prod(u)^((R.θ) / (1 - R.θ))
-
-    logpdf_value = log(term_numerator) - log(term_denominator) + log(term_product)
-
-    return logpdf_value
+    l_den = (d-1) * log(1-R.θ) + log(d + R.θ -1)
+    l_num = log(d - 1 + R.θ * u_ordered[d]^((1 - R.θ - d) / (1 - R.θ)))
+    l_prd = (R.θ) / (1 - R.θ) * log(prod(u))
+    return l_num - l_den + l_prd
 end
-
 function Distributions._rand!(rng::Distributions.AbstractRNG, R::RafteryCopula{d,P}, x::AbstractVector{T}) where {d,P,T <: Real}
-
-    dim = length(x)
-
     # Step 1: Generate independent values u, u_1, ..., u_d from a uniform distribution [0, 1]
-    u = rand(rng, dim+1)
-    
+    u = rand(rng, d)
+
     # Step 2: Generate j from a Bernoulli distribution with parameter θ
-    j = rand(Distributions.Bernoulli(R.θ), 1)
-    uj = u[1]^j[1]
+    uj = (rand(rng)<R.θ) ? rand(rng) : 1
+
     # Step 3: Calculate v_1, ..., v_d
-    for i in 2:dim+1
-        x[i-1] = u[i]^(1 - R.θ) * uj
+    for i in 1:d
+        x[i] = u[i]^(1 - R.θ) * uj
     end
 
     return x

test/RafteryTest.jl

@@ -13,44 +13,23 @@ end
     rng = StableRNG(123)
     for d in [2, 3, 4]
         F = RafteryCopula(d, 0.5)
-
-        # Test CDF with some random values
-        u = rand(d)
-        cdf_value = cdf(F, u)
+        cdf_value = cdf(F, rand(d))
+        pdf_value = pdf(F,rand(d))
         @test cdf_value >= 0 && cdf_value <= 1
+        @test pdf_value >= 0 
     end
 
     @test cdf(RafteryCopula(2, 0.8), [0.2, 0.5]) ≈ 0.199432 atol=1e-5
     @test cdf(RafteryCopula(2, 0.5), [0.3, 0.8]) ≈ 0.2817 atol=1e-5
-    @test_broken cdf(RafteryCopula(3, 0.5), [0.1, 0.2, 0.3]) ≈ 0.08325884 atol=1e-5
-    @test_broken cdf(RafteryCopula(3, 0.1), [0.4, 0.8, 0.2]) ≈ 0.120415 atol=1e-5    
-end
-
-@testitem "RafteryCopula PDF" begin
-    using StableRNGs
-    using Distributions
-    rng = StableRNG(123)
-    for d in [2, 3, 4]
-        F = RafteryCopula(d, 0.5)
+    @test cdf(RafteryCopula(3, 0.5), [0.1, 0.2, 0.3]) ≈ 0.08236007 atol=1e-5
+    @test cdf(RafteryCopula(3, 0.1), [0.4, 0.8, 0.2]) ≈ 0.08581997 atol=1e-5    
 
-        # Test PDF with some random values
-        u = rand(rng,d)
-        @test_broken 0 <= pdf(F, u) <= 1
-    end
-    examples = [
-        ([0.2, 0.5], [0.114055555, 1e-4], 0.8),
-        ([0.3, 0.8], [0.6325, 1e-4], 0.5),
-        ([0.1, 0.2, 0.3], [1.9945086, 1e-4], 0.5),
-        ([0.4, 0.8, 0.2], [0.939229, 1e-4], 0.1),
-    ]
-
-    for (u, expected, θ) in examples
-        copula = RafteryCopula(length(u), θ)
-        @test_broken pdf(copula, u) ≈ expected[1] atol=expected[2]
-    end
+    @test pdf(RafteryCopula(2, 0.8), [0.2, 0.5]) ≈ 0.114055555 atol=1e-4
+    @test pdf(RafteryCopula(2, 0.5), [0.3, 0.8]) ≈ 0.6325 atol=1e-4
+    @test pdf(RafteryCopula(3, 0.5), [0.1, 0.2, 0.3]) ≈ 1.9945086 atol=1e-4
+    @test pdf(RafteryCopula(3, 0.1), [0.4, 0.8, 0.2]) ≈ 0.939229 atol=1e-4
 end
 
-
 @testitem "RafteryCopula Sampling" begin
     using StableRNGs
     rng = StableRNG(123)
@@ -59,3 +38,71 @@ end
     samples = rand(rng,F, n_samples)
     @test size(samples) == (3, n_samples)
 end
+
+@testitem "Check against manual version - CDF" begin
+    # https://github.com/lrnv/Copulas.jl/pull/137
+    using Distributions
+    function prueba_CDF(R::Vector{T}, u::Vector{T}) where T
+        # Order the vector u
+        θ = R[1]
+        # println("param:", θ)
+        d = round(Int, R[2])
+        # println("dimension:", d)
+        u_ordered = sort(u)
+        # println("Sorted vector: ", u_ordered)
+
+        term1 = u_ordered[1]
+        # println("Term 1: ", term1)
+
+        term2 = ((1 - θ) * (1 -d)) / (1 - θ - d) * prod(u)^(1/(1 - θ))
+        # println("Term 2: ", term2)
+
+        term3 = 0.0
+        for i in 2:d # <<<<<<<<--------------- This should be 2:d and not 2:length(R) since length(R) is not the dimension. 
+            prod_prev = prod(u_ordered[1:i-1])
+            term3_part = ((θ * (1 - θ)) / ((1 - θ - i) * (2 - θ - i))) * prod_prev^(1/(1 - θ)) * u_ordered[i]^((2 - θ - i) / (1 - θ))
+            # println("Term 3 (part $i): ", term3_part)
+            term3 += term3_part
+        end
+
+        # Combine the terms to get the cumulative distribution function
+        cdf_value = term1 + term2 - term3
+        # println("Final CDF value: ", cdf_value)
+
+        return cdf_value
+    end
+    @test prueba_CDF([0.5,3], [0.1,0.2,0.3]) ≈ 0.08236 atol=1e-4 # According to https://github.com/lrnv/Copulas.jl/pull/137#issuecomment-1953365273
+    @test prueba_CDF([0.5,3], [0.1,0.2,0.3]) ≈ cdf(RafteryCopula(3,0.5), [0.1,0.2,0.3])
+    @test prueba_CDF([0.8,2], [0.1,0.2]) ≈ cdf(RafteryCopula(2,0.8), [0.1,0.2])
+    @test prueba_CDF([0.2,2], [0.8,0.2]) ≈ cdf(RafteryCopula(2,0.2), [0.8,0.2])
+end
+
+@testitem "Check against manual version - PDF" begin
+    # https://github.com/lrnv/Copulas.jl/pull/137
+    using Distributions
+    function prueba_PDF(R::Vector{T}, u::Vector{T}) where T
+        # Order the vector u
+        θ = R[1]
+        d = round(Int, R[2])
+        u_ordered = sort(u)
+        # println("Sorted vector: ", u_ordered)
+
+        term1 = (1/(((1-θ)^(d-1))*(1-θ-d)))
+        # println("Term 1: ", term1)
+
+        term2 = (1-d-θ*(u_ordered[d])^((1-θ-d)/(1-θ)))
+        # println("Term 2: ", term2)
+
+        term3 = (prod(u))^((θ)/(1-θ))
+        # println(term3)
+        # Combine the terms to get the density distribution function
+        pdf_value = term1*term2*term3
+        # println("Final PDF value: ", pdf_value)
+
+        return pdf_value
+    end
+    @test prueba_PDF([0.5,3], [0.1,0.2,0.3]) ≈ 1.99450 atol=1e-4 # According to https://github.com/lrnv/Copulas.jl/pull/137#issuecomment-1953375141
+    @test prueba_PDF([0.5,3], [0.1,0.2,0.3]) ≈ pdf(RafteryCopula(3,0.5), [0.1,0.2,0.3])
+    @test prueba_PDF([0.8,2], [0.1,0.2]) ≈ pdf(RafteryCopula(2,0.8), [0.1,0.2])
+    @test prueba_PDF([0.2,2], [0.8,0.2]) ≈ pdf(RafteryCopula(2,0.2), [0.8,0.2])
+end

test/kendall_tau_notnan.jl

@@ -32,7 +32,6 @@
     )
 
     for C in cops
-        @show C
         @test !isnan(Copulas.τ(C))
     end
     @test_broken Copulas.τ(ArchimedeanCopula(2,i𝒲(LogNormal(),2))) # not implemented. 

test/margins_uniformity.jl

@@ -94,7 +94,6 @@
     n = 1000
     U = Uniform(0,1)
     for C in cops
-        @show C
         d = length(C)
         CT = typeof(C)
         rng = StableRNG(123)