From e227918496b91f53dfed3491d2904925cbac90c5 Mon Sep 17 00:00:00 2001 From: Santiago Jimenez Ramos Date: Thu, 26 Oct 2023 00:56:37 -0300 Subject: [PATCH 1/7] The file PlackettCopula.jl is added that includes: an instance of the bivariate Plackett Copula, calculation of the CDF and PDF functions, the random sample generator and the Spearman Rho calculation. Additionally, I added some testing these functions in biv_plackett.jl. --- src/Copulas.jl | 4 +- src/MiscellaneousCopulas/PlackettCopula.jl | 95 ++++++++++++++++++++++ test/biv_Plackett.jl | 40 +++++++++ 3 files changed, 138 insertions(+), 1 deletion(-) create mode 100644 src/MiscellaneousCopulas/PlackettCopula.jl create mode 100644 test/biv_Plackett.jl diff --git a/src/Copulas.jl b/src/Copulas.jl index 6b16b457..6736174e 100644 --- a/src/Copulas.jl +++ b/src/Copulas.jl @@ -60,9 +60,11 @@ end Copulas include("MiscellaneousCopulas/MCopula.jl") include("MiscellaneousCopulas/WCopula.jl") include("MiscellaneousCopulas/SurvivalCopula.jl") + include("MiscellaneousCopulas/PlackettCopula.jl") export MCopula, Wcopula, - SurvivalCopula + SurvivalCopula, + PlackettCopula # PrecompileTools stuff @setup_workload begin diff --git a/src/MiscellaneousCopulas/PlackettCopula.jl b/src/MiscellaneousCopulas/PlackettCopula.jl new file mode 100644 index 00000000..18c70865 --- /dev/null +++ b/src/MiscellaneousCopulas/PlackettCopula.jl @@ -0,0 +1,95 @@ +using Distributions, Copulas +#= Details about Plackett copulation are found in Joe, H. (2014). + Dependence modeling with copulas. CRC press, Page.164 +==# + +#Create an instance of the Plackett copula +struct PlackettCopula{P} <: ContinuousMultivariateDistribution + θ::P # Copula parameter + + function PlackettCopula(θ) where {P} + if θ == 1.0 + return IndependentCopula(2) + elseif θ == 0 + return MCopula(2) + elseif θ == Inf + return WCopula(2) + else + θ >= 0 || throw(ArgumentError("Theta must be non-negative")) + return new{typeof(θ)}(θ) + end + end +end + +Base.length(S::PlackettCopula{P}) where {P} = 1 +Base.eltype(S::PlackettCopula{P}) where {P} = Float64 + +import Distributions: cdf, pdf + +# CDF calculation for bivariate Plackett Copula +function cdf(S::PlackettCopula{P}, uv::Tuple{Float64, Float64}) where {P} + u, v = uv + η = S.θ - 1 + + if S.θ == 1.0 + return cdf(IndependentCopula(2), uv) + elseif S.θ == Inf + return cdf(WCopula(2), uv) + elseif S.θ == 0 + return cdf(MCopula(2), uv) + else + term1 = 1 + η * (u + v) + term2 = sqrt(term1^2 - 4 * S.θ * η * u * v) + return 0.5 * η^(-1) * (term1 - term2) + end +end + +# PDF calculation for bivariate Plackett Copula +function pdf(S::PlackettCopula{P}, uv::Tuple{Float64, Float64}) where {P} + u, v = uv + η = S.θ - 1 + + if S.θ == 1.0 + return pdf(IndependentCopula(2), uv) + elseif S.θ == Inf + throw(ArgumentError("Indefinite density")) + elseif S.θ == 0 + throw(ArgumentError("Indefinite density")) + else + term1 = S.θ * (1 + η * (u + v - 2 * u * v)) + term2 = (1+η*(u+v))^2-4*(S.θ)*η*u*v + return term1/term2^(3/2) + end +end +import Random + +#= Details about the algorithm to generate copula samples + can be seen in the following references + Johnson, Mark E. Multivariate statistical simulation: + A guide to selecting and generating continuous multivariate distributions. + Vol. 192. John Wiley & Sons, 1987. Page 193. + Nelsen, Roger B. An introduction to copulas. Springer, 2006. Exercise 3.38. +==# + +# Copula random sample simulator +function rand(rng::Random.AbstractRNG, c::PlackettCopula{P}, n::Int) where P + samples = Matrix{Float64}(undef, 2, n) + for i in 1:n + u = rand(rng) + t = rand(rng) + a = t * (1 - t) + b = c.θ + a * (c.θ - 1)^2 + cc = 2a * (u * c.θ^2 + 1 - u) + c.θ * (1 - 2a) + d = sqrt(c.θ) * sqrt(c.θ + 4a * u * (1 - u) * (1 - c.θ)^2) + v = (cc - (1 - 2t) * d) / (2b) + samples[1, i] = u + samples[2, i] = v + end + return samples +end + +# Calculate Spearman's rho based on the PlackettCopula parameters +function spearman_rho(c::PlackettCopula{P}) where P + rho = (c.θ+1)/(c.θ-1)-(2*c.θ*log(c.θ)/(c.θ-1)^2) + return rho +end \ No newline at end of file diff --git a/test/biv_Plackett.jl b/test/biv_Plackett.jl new file mode 100644 index 00000000..1168dc0a --- /dev/null +++ b/test/biv_Plackett.jl @@ -0,0 +1,40 @@ +@testset "PlackettCopula Tests" begin + @testset "Constructor" begin + @test_broken isa(PlackettCopula(1), Independence) + @test_broken isa(PlackettCopula(Inf),WCopula) # should work in any dimenisons if theta is smaller than the bound. + @test_broken isa(PlackettCopula(0),MCopula) + @test_throws ArgumentError PlackettCopula(-0.5) + end + + @testset "CDF" begin + u = 0.1:0.18:1 + v = 0.4:0.1:0.9 + l1 = [0.0553778 ,0.1743884 ,0.3166277 ,0.4823228 ,0.6743114 ,0.9 ] + l2 = [0.02620873, 0.1056116, 0.2349113, 0.4162573, 0.6419255 ,0.9] + + for i in 1:6 + @test cdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i] + @test cdf(PlackettCopula(0.5), [u[i], v[i]]) ≈ l2[i] + end + end + + @testset "PDF" begin + u = 0.1:0.18:1 + v = 0.4:0.1:0.9 + l1 = [1.059211 ,1.023291 ,1.038467 ,1.110077 ,1.272959 ,1.652893 ] + l2 = [0.8446203 ,1.023291, 1.064891, 0.9360171, 0.7346612, 0.5540166] + + for i in 1:6 + @test pdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i] + @test pdf(PlackettCopula(0.5), [u[i], v[i]]) ≈ l2[i] + end + end + + @testset "Sampling" begin + using Random + n_samples = 100 + c = PlackettCopula(0.8) + samples = rand(c, n_samples) + @test size(samples) == (2, n_samples) + end +end From 0a51c443c72f9834959fc16b26e5273705067c9e Mon Sep 17 00:00:00 2001 From: Oskar Laverny Date: Thu, 26 Oct 2023 14:56:14 +0200 Subject: [PATCH 2/7] first shot --- src/MiscellaneousCopulas/PlackettCopula.jl | 89 +++++++++------------- test/biv_Plackett.jl | 62 ++++++++------- 2 files changed, 68 insertions(+), 83 deletions(-) diff --git a/src/MiscellaneousCopulas/PlackettCopula.jl b/src/MiscellaneousCopulas/PlackettCopula.jl index 18c70865..f174436c 100644 --- a/src/MiscellaneousCopulas/PlackettCopula.jl +++ b/src/MiscellaneousCopulas/PlackettCopula.jl @@ -1,14 +1,13 @@ -using Distributions, Copulas #= Details about Plackett copulation are found in Joe, H. (2014). Dependence modeling with copulas. CRC press, Page.164 ==# #Create an instance of the Plackett copula -struct PlackettCopula{P} <: ContinuousMultivariateDistribution +struct PlackettCopula{P} <: Copula{2} # since it is only bivariate. θ::P # Copula parameter function PlackettCopula(θ) where {P} - if θ == 1.0 + if θ == 1 return IndependentCopula(2) elseif θ == 0 return MCopula(2) @@ -21,45 +20,25 @@ struct PlackettCopula{P} <: ContinuousMultivariateDistribution end end -Base.length(S::PlackettCopula{P}) where {P} = 1 -Base.eltype(S::PlackettCopula{P}) where {P} = Float64 - -import Distributions: cdf, pdf +# Base.length(S::PlackettCopula{P}) where {P} = 2 # length should return the dimension of the copula, bnut i think it is already working without this definition. +Base.eltype(S::PlackettCopula{P}) where {P} = P # this shuold be P. # CDF calculation for bivariate Plackett Copula -function cdf(S::PlackettCopula{P}, uv::Tuple{Float64, Float64}) where {P} +function Distributions.cdf(S::PlackettCopula{P}, uv) where {P} u, v = uv η = S.θ - 1 - - if S.θ == 1.0 - return cdf(IndependentCopula(2), uv) - elseif S.θ == Inf - return cdf(WCopula(2), uv) - elseif S.θ == 0 - return cdf(MCopula(2), uv) - else - term1 = 1 + η * (u + v) - term2 = sqrt(term1^2 - 4 * S.θ * η * u * v) - return 0.5 * η^(-1) * (term1 - term2) - end + term1 = 1 + η * (u + v) + term2 = sqrt(term1^2 - 4 * S.θ * η * u * v) + return 0.5 * η^(-1) * (term1 - term2) end # PDF calculation for bivariate Plackett Copula -function pdf(S::PlackettCopula{P}, uv::Tuple{Float64, Float64}) where {P} +function Distributions._logpdf(S::PlackettCopula{P}, uv) where {P} u, v = uv η = S.θ - 1 - - if S.θ == 1.0 - return pdf(IndependentCopula(2), uv) - elseif S.θ == Inf - throw(ArgumentError("Indefinite density")) - elseif S.θ == 0 - throw(ArgumentError("Indefinite density")) - else - term1 = S.θ * (1 + η * (u + v - 2 * u * v)) - term2 = (1+η*(u+v))^2-4*(S.θ)*η*u*v - return term1/term2^(3/2) - end + term1 = S.θ * (1 + η * (u + v - 2 * u * v)) + term2 = (1+η*(u+v))^2-4*(S.θ)*η*u*v + return log(term1) - 3 * log(term2)/2 # since we are supposed to return the logpdf. end import Random @@ -71,25 +50,33 @@ import Random Nelsen, Roger B. An introduction to copulas. Springer, 2006. Exercise 3.38. ==# -# Copula random sample simulator -function rand(rng::Random.AbstractRNG, c::PlackettCopula{P}, n::Int) where P - samples = Matrix{Float64}(undef, 2, n) - for i in 1:n - u = rand(rng) - t = rand(rng) - a = t * (1 - t) - b = c.θ + a * (c.θ - 1)^2 - cc = 2a * (u * c.θ^2 + 1 - u) + c.θ * (1 - 2a) - d = sqrt(c.θ) * sqrt(c.θ + 4a * u * (1 - u) * (1 - c.θ)^2) - v = (cc - (1 - 2t) * d) / (2b) - samples[1, i] = u - samples[2, i] = v - end - return samples +function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, x::AbstractVector{T}) where {T<:Real, CT<:PlackettCopula} + u = rand(rng) + t = rand(rng) + a = t * (1 - t) + b = c.θ + a * (c.θ - 1)^2 + cc = 2a * (u * c.θ^2 + 1 - u) + c.θ * (1 - 2a) + d = sqrt(c.θ) * sqrt(c.θ + 4a * u * (1 - u) * (1 - c.θ)^2) + v = (cc - (1 - 2t) * d) / (2b) + x[1] = u + x[2] = v + return x +end +function Base.rand(rng::Distributions.AbstractRNG,C::CT) where CT<: PlackettCopula + x = rand(rng,length(C)) + u = rand(rng) + t = rand(rng) + a = t * (1 - t) + b = c.θ + a * (c.θ - 1)^2 + cc = 2a * (u * c.θ^2 + 1 - u) + c.θ * (1 - 2a) + d = sqrt(c.θ) * sqrt(c.θ + 4a * u * (1 - u) * (1 - c.θ)^2) + v = (cc - (1 - 2t) * d) / (2b) + x[1] = u + x[2] = v + return x end # Calculate Spearman's rho based on the PlackettCopula parameters -function spearman_rho(c::PlackettCopula{P}) where P - rho = (c.θ+1)/(c.θ-1)-(2*c.θ*log(c.θ)/(c.θ-1)^2) - return rho +function ρ(c::PlackettCopula{P}) where P + return (c.θ+1)/(c.θ-1)-(2*c.θ*log(c.θ)/(c.θ-1)^2) end \ No newline at end of file diff --git a/test/biv_Plackett.jl b/test/biv_Plackett.jl index 1168dc0a..c03dc960 100644 --- a/test/biv_Plackett.jl +++ b/test/biv_Plackett.jl @@ -1,40 +1,38 @@ -@testset "PlackettCopula Tests" begin - @testset "Constructor" begin - @test_broken isa(PlackettCopula(1), Independence) - @test_broken isa(PlackettCopula(Inf),WCopula) # should work in any dimenisons if theta is smaller than the bound. - @test_broken isa(PlackettCopula(0),MCopula) - @test_throws ArgumentError PlackettCopula(-0.5) - end +@testitem "PlackettCopula Constructor" begin + @test_broken isa(PlackettCopula(1), Independence) + @test_broken isa(PlackettCopula(Inf),WCopula) # should work in any dimenisons if theta is smaller than the bound. + @test_broken isa(PlackettCopula(0),MCopula) + @test_throws ArgumentError PlackettCopula(-0.5) +end - @testset "CDF" begin - u = 0.1:0.18:1 - v = 0.4:0.1:0.9 - l1 = [0.0553778 ,0.1743884 ,0.3166277 ,0.4823228 ,0.6743114 ,0.9 ] - l2 = [0.02620873, 0.1056116, 0.2349113, 0.4162573, 0.6419255 ,0.9] +@testitem "PlackettCopula CDF" begin + u = 0.1:0.18:1 + v = 0.4:0.1:0.9 + l1 = [0.0553778 ,0.1743884 ,0.3166277 ,0.4823228 ,0.6743114 ,0.9 ] + l2 = [0.02620873, 0.1056116, 0.2349113, 0.4162573, 0.6419255 ,0.9] - for i in 1:6 - @test cdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i] - @test cdf(PlackettCopula(0.5), [u[i], v[i]]) ≈ l2[i] - end + for i in 1:6 + @test cdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i] + @test cdf(PlackettCopula(0.5), [u[i], v[i]]) ≈ l2[i] end +end - @testset "PDF" begin - u = 0.1:0.18:1 - v = 0.4:0.1:0.9 - l1 = [1.059211 ,1.023291 ,1.038467 ,1.110077 ,1.272959 ,1.652893 ] - l2 = [0.8446203 ,1.023291, 1.064891, 0.9360171, 0.7346612, 0.5540166] +@testitem "PlackettCopula PDF" begin + u = 0.1:0.18:1 + v = 0.4:0.1:0.9 + l1 = [1.059211 ,1.023291 ,1.038467 ,1.110077 ,1.272959 ,1.652893 ] + l2 = [0.8446203 ,1.023291, 1.064891, 0.9360171, 0.7346612, 0.5540166] - for i in 1:6 - @test pdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i] - @test pdf(PlackettCopula(0.5), [u[i], v[i]]) ≈ l2[i] - end + for i in 1:6 + @test pdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i] + @test pdf(PlackettCopula(0.5), [u[i], v[i]]) ≈ l2[i] end +end - @testset "Sampling" begin - using Random - n_samples = 100 - c = PlackettCopula(0.8) - samples = rand(c, n_samples) - @test size(samples) == (2, n_samples) - end +@testitem "PlackettCopula Sampling" begin + using Random + n_samples = 100 + c = PlackettCopula(0.8) + samples = rand(c, n_samples) + @test size(samples) == (2, n_samples) end From 77f8e8fd2fbb3a1d122f8d2a403da340025e28ee Mon Sep 17 00:00:00 2001 From: Santiago Jimenez Ramos Date: Thu, 26 Oct 2023 00:56:37 -0300 Subject: [PATCH 3/7] The file PlackettCopula.jl is added that includes: an instance of the bivariate Plackett Copula, calculation of the CDF and PDF functions, the random sample generator and the Spearman Rho calculation. Additionally, I added some testing these functions in biv_plackett.jl. --- src/Copulas.jl | 3 +- src/MiscellaneousCopulas/PlackettCopula.jl | 95 ++++++++++++++++++++++ test/biv_Plackett.jl | 40 +++++++++ 3 files changed, 137 insertions(+), 1 deletion(-) create mode 100644 src/MiscellaneousCopulas/PlackettCopula.jl create mode 100644 test/biv_Plackett.jl diff --git a/src/Copulas.jl b/src/Copulas.jl index 2f25dbf3..fa2f6e79 100644 --- a/src/Copulas.jl +++ b/src/Copulas.jl @@ -59,8 +59,9 @@ end Copulas include("MiscellaneousCopulas/MCopula.jl") include("MiscellaneousCopulas/WCopula.jl") include("MiscellaneousCopulas/SurvivalCopula.jl") + include("MiscellaneousCopulas/PlackettCopula.jl") export MCopula, WCopula, SurvivalCopula - + PlackettCopula end diff --git a/src/MiscellaneousCopulas/PlackettCopula.jl b/src/MiscellaneousCopulas/PlackettCopula.jl new file mode 100644 index 00000000..18c70865 --- /dev/null +++ b/src/MiscellaneousCopulas/PlackettCopula.jl @@ -0,0 +1,95 @@ +using Distributions, Copulas +#= Details about Plackett copulation are found in Joe, H. (2014). + Dependence modeling with copulas. CRC press, Page.164 +==# + +#Create an instance of the Plackett copula +struct PlackettCopula{P} <: ContinuousMultivariateDistribution + θ::P # Copula parameter + + function PlackettCopula(θ) where {P} + if θ == 1.0 + return IndependentCopula(2) + elseif θ == 0 + return MCopula(2) + elseif θ == Inf + return WCopula(2) + else + θ >= 0 || throw(ArgumentError("Theta must be non-negative")) + return new{typeof(θ)}(θ) + end + end +end + +Base.length(S::PlackettCopula{P}) where {P} = 1 +Base.eltype(S::PlackettCopula{P}) where {P} = Float64 + +import Distributions: cdf, pdf + +# CDF calculation for bivariate Plackett Copula +function cdf(S::PlackettCopula{P}, uv::Tuple{Float64, Float64}) where {P} + u, v = uv + η = S.θ - 1 + + if S.θ == 1.0 + return cdf(IndependentCopula(2), uv) + elseif S.θ == Inf + return cdf(WCopula(2), uv) + elseif S.θ == 0 + return cdf(MCopula(2), uv) + else + term1 = 1 + η * (u + v) + term2 = sqrt(term1^2 - 4 * S.θ * η * u * v) + return 0.5 * η^(-1) * (term1 - term2) + end +end + +# PDF calculation for bivariate Plackett Copula +function pdf(S::PlackettCopula{P}, uv::Tuple{Float64, Float64}) where {P} + u, v = uv + η = S.θ - 1 + + if S.θ == 1.0 + return pdf(IndependentCopula(2), uv) + elseif S.θ == Inf + throw(ArgumentError("Indefinite density")) + elseif S.θ == 0 + throw(ArgumentError("Indefinite density")) + else + term1 = S.θ * (1 + η * (u + v - 2 * u * v)) + term2 = (1+η*(u+v))^2-4*(S.θ)*η*u*v + return term1/term2^(3/2) + end +end +import Random + +#= Details about the algorithm to generate copula samples + can be seen in the following references + Johnson, Mark E. Multivariate statistical simulation: + A guide to selecting and generating continuous multivariate distributions. + Vol. 192. John Wiley & Sons, 1987. Page 193. + Nelsen, Roger B. An introduction to copulas. Springer, 2006. Exercise 3.38. +==# + +# Copula random sample simulator +function rand(rng::Random.AbstractRNG, c::PlackettCopula{P}, n::Int) where P + samples = Matrix{Float64}(undef, 2, n) + for i in 1:n + u = rand(rng) + t = rand(rng) + a = t * (1 - t) + b = c.θ + a * (c.θ - 1)^2 + cc = 2a * (u * c.θ^2 + 1 - u) + c.θ * (1 - 2a) + d = sqrt(c.θ) * sqrt(c.θ + 4a * u * (1 - u) * (1 - c.θ)^2) + v = (cc - (1 - 2t) * d) / (2b) + samples[1, i] = u + samples[2, i] = v + end + return samples +end + +# Calculate Spearman's rho based on the PlackettCopula parameters +function spearman_rho(c::PlackettCopula{P}) where P + rho = (c.θ+1)/(c.θ-1)-(2*c.θ*log(c.θ)/(c.θ-1)^2) + return rho +end \ No newline at end of file diff --git a/test/biv_Plackett.jl b/test/biv_Plackett.jl new file mode 100644 index 00000000..1168dc0a --- /dev/null +++ b/test/biv_Plackett.jl @@ -0,0 +1,40 @@ +@testset "PlackettCopula Tests" begin + @testset "Constructor" begin + @test_broken isa(PlackettCopula(1), Independence) + @test_broken isa(PlackettCopula(Inf),WCopula) # should work in any dimenisons if theta is smaller than the bound. + @test_broken isa(PlackettCopula(0),MCopula) + @test_throws ArgumentError PlackettCopula(-0.5) + end + + @testset "CDF" begin + u = 0.1:0.18:1 + v = 0.4:0.1:0.9 + l1 = [0.0553778 ,0.1743884 ,0.3166277 ,0.4823228 ,0.6743114 ,0.9 ] + l2 = [0.02620873, 0.1056116, 0.2349113, 0.4162573, 0.6419255 ,0.9] + + for i in 1:6 + @test cdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i] + @test cdf(PlackettCopula(0.5), [u[i], v[i]]) ≈ l2[i] + end + end + + @testset "PDF" begin + u = 0.1:0.18:1 + v = 0.4:0.1:0.9 + l1 = [1.059211 ,1.023291 ,1.038467 ,1.110077 ,1.272959 ,1.652893 ] + l2 = [0.8446203 ,1.023291, 1.064891, 0.9360171, 0.7346612, 0.5540166] + + for i in 1:6 + @test pdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i] + @test pdf(PlackettCopula(0.5), [u[i], v[i]]) ≈ l2[i] + end + end + + @testset "Sampling" begin + using Random + n_samples = 100 + c = PlackettCopula(0.8) + samples = rand(c, n_samples) + @test size(samples) == (2, n_samples) + end +end From 9b126c3b1405704d20da2d9bd5d3427fef5979ad Mon Sep 17 00:00:00 2001 From: Oskar Laverny Date: Thu, 26 Oct 2023 14:56:14 +0200 Subject: [PATCH 4/7] first shot --- src/MiscellaneousCopulas/PlackettCopula.jl | 89 +++++++++------------- test/biv_Plackett.jl | 62 ++++++++------- 2 files changed, 68 insertions(+), 83 deletions(-) diff --git a/src/MiscellaneousCopulas/PlackettCopula.jl b/src/MiscellaneousCopulas/PlackettCopula.jl index 18c70865..f174436c 100644 --- a/src/MiscellaneousCopulas/PlackettCopula.jl +++ b/src/MiscellaneousCopulas/PlackettCopula.jl @@ -1,14 +1,13 @@ -using Distributions, Copulas #= Details about Plackett copulation are found in Joe, H. (2014). Dependence modeling with copulas. CRC press, Page.164 ==# #Create an instance of the Plackett copula -struct PlackettCopula{P} <: ContinuousMultivariateDistribution +struct PlackettCopula{P} <: Copula{2} # since it is only bivariate. θ::P # Copula parameter function PlackettCopula(θ) where {P} - if θ == 1.0 + if θ == 1 return IndependentCopula(2) elseif θ == 0 return MCopula(2) @@ -21,45 +20,25 @@ struct PlackettCopula{P} <: ContinuousMultivariateDistribution end end -Base.length(S::PlackettCopula{P}) where {P} = 1 -Base.eltype(S::PlackettCopula{P}) where {P} = Float64 - -import Distributions: cdf, pdf +# Base.length(S::PlackettCopula{P}) where {P} = 2 # length should return the dimension of the copula, bnut i think it is already working without this definition. +Base.eltype(S::PlackettCopula{P}) where {P} = P # this shuold be P. # CDF calculation for bivariate Plackett Copula -function cdf(S::PlackettCopula{P}, uv::Tuple{Float64, Float64}) where {P} +function Distributions.cdf(S::PlackettCopula{P}, uv) where {P} u, v = uv η = S.θ - 1 - - if S.θ == 1.0 - return cdf(IndependentCopula(2), uv) - elseif S.θ == Inf - return cdf(WCopula(2), uv) - elseif S.θ == 0 - return cdf(MCopula(2), uv) - else - term1 = 1 + η * (u + v) - term2 = sqrt(term1^2 - 4 * S.θ * η * u * v) - return 0.5 * η^(-1) * (term1 - term2) - end + term1 = 1 + η * (u + v) + term2 = sqrt(term1^2 - 4 * S.θ * η * u * v) + return 0.5 * η^(-1) * (term1 - term2) end # PDF calculation for bivariate Plackett Copula -function pdf(S::PlackettCopula{P}, uv::Tuple{Float64, Float64}) where {P} +function Distributions._logpdf(S::PlackettCopula{P}, uv) where {P} u, v = uv η = S.θ - 1 - - if S.θ == 1.0 - return pdf(IndependentCopula(2), uv) - elseif S.θ == Inf - throw(ArgumentError("Indefinite density")) - elseif S.θ == 0 - throw(ArgumentError("Indefinite density")) - else - term1 = S.θ * (1 + η * (u + v - 2 * u * v)) - term2 = (1+η*(u+v))^2-4*(S.θ)*η*u*v - return term1/term2^(3/2) - end + term1 = S.θ * (1 + η * (u + v - 2 * u * v)) + term2 = (1+η*(u+v))^2-4*(S.θ)*η*u*v + return log(term1) - 3 * log(term2)/2 # since we are supposed to return the logpdf. end import Random @@ -71,25 +50,33 @@ import Random Nelsen, Roger B. An introduction to copulas. Springer, 2006. Exercise 3.38. ==# -# Copula random sample simulator -function rand(rng::Random.AbstractRNG, c::PlackettCopula{P}, n::Int) where P - samples = Matrix{Float64}(undef, 2, n) - for i in 1:n - u = rand(rng) - t = rand(rng) - a = t * (1 - t) - b = c.θ + a * (c.θ - 1)^2 - cc = 2a * (u * c.θ^2 + 1 - u) + c.θ * (1 - 2a) - d = sqrt(c.θ) * sqrt(c.θ + 4a * u * (1 - u) * (1 - c.θ)^2) - v = (cc - (1 - 2t) * d) / (2b) - samples[1, i] = u - samples[2, i] = v - end - return samples +function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, x::AbstractVector{T}) where {T<:Real, CT<:PlackettCopula} + u = rand(rng) + t = rand(rng) + a = t * (1 - t) + b = c.θ + a * (c.θ - 1)^2 + cc = 2a * (u * c.θ^2 + 1 - u) + c.θ * (1 - 2a) + d = sqrt(c.θ) * sqrt(c.θ + 4a * u * (1 - u) * (1 - c.θ)^2) + v = (cc - (1 - 2t) * d) / (2b) + x[1] = u + x[2] = v + return x +end +function Base.rand(rng::Distributions.AbstractRNG,C::CT) where CT<: PlackettCopula + x = rand(rng,length(C)) + u = rand(rng) + t = rand(rng) + a = t * (1 - t) + b = c.θ + a * (c.θ - 1)^2 + cc = 2a * (u * c.θ^2 + 1 - u) + c.θ * (1 - 2a) + d = sqrt(c.θ) * sqrt(c.θ + 4a * u * (1 - u) * (1 - c.θ)^2) + v = (cc - (1 - 2t) * d) / (2b) + x[1] = u + x[2] = v + return x end # Calculate Spearman's rho based on the PlackettCopula parameters -function spearman_rho(c::PlackettCopula{P}) where P - rho = (c.θ+1)/(c.θ-1)-(2*c.θ*log(c.θ)/(c.θ-1)^2) - return rho +function ρ(c::PlackettCopula{P}) where P + return (c.θ+1)/(c.θ-1)-(2*c.θ*log(c.θ)/(c.θ-1)^2) end \ No newline at end of file diff --git a/test/biv_Plackett.jl b/test/biv_Plackett.jl index 1168dc0a..c03dc960 100644 --- a/test/biv_Plackett.jl +++ b/test/biv_Plackett.jl @@ -1,40 +1,38 @@ -@testset "PlackettCopula Tests" begin - @testset "Constructor" begin - @test_broken isa(PlackettCopula(1), Independence) - @test_broken isa(PlackettCopula(Inf),WCopula) # should work in any dimenisons if theta is smaller than the bound. - @test_broken isa(PlackettCopula(0),MCopula) - @test_throws ArgumentError PlackettCopula(-0.5) - end +@testitem "PlackettCopula Constructor" begin + @test_broken isa(PlackettCopula(1), Independence) + @test_broken isa(PlackettCopula(Inf),WCopula) # should work in any dimenisons if theta is smaller than the bound. + @test_broken isa(PlackettCopula(0),MCopula) + @test_throws ArgumentError PlackettCopula(-0.5) +end - @testset "CDF" begin - u = 0.1:0.18:1 - v = 0.4:0.1:0.9 - l1 = [0.0553778 ,0.1743884 ,0.3166277 ,0.4823228 ,0.6743114 ,0.9 ] - l2 = [0.02620873, 0.1056116, 0.2349113, 0.4162573, 0.6419255 ,0.9] +@testitem "PlackettCopula CDF" begin + u = 0.1:0.18:1 + v = 0.4:0.1:0.9 + l1 = [0.0553778 ,0.1743884 ,0.3166277 ,0.4823228 ,0.6743114 ,0.9 ] + l2 = [0.02620873, 0.1056116, 0.2349113, 0.4162573, 0.6419255 ,0.9] - for i in 1:6 - @test cdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i] - @test cdf(PlackettCopula(0.5), [u[i], v[i]]) ≈ l2[i] - end + for i in 1:6 + @test cdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i] + @test cdf(PlackettCopula(0.5), [u[i], v[i]]) ≈ l2[i] end +end - @testset "PDF" begin - u = 0.1:0.18:1 - v = 0.4:0.1:0.9 - l1 = [1.059211 ,1.023291 ,1.038467 ,1.110077 ,1.272959 ,1.652893 ] - l2 = [0.8446203 ,1.023291, 1.064891, 0.9360171, 0.7346612, 0.5540166] +@testitem "PlackettCopula PDF" begin + u = 0.1:0.18:1 + v = 0.4:0.1:0.9 + l1 = [1.059211 ,1.023291 ,1.038467 ,1.110077 ,1.272959 ,1.652893 ] + l2 = [0.8446203 ,1.023291, 1.064891, 0.9360171, 0.7346612, 0.5540166] - for i in 1:6 - @test pdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i] - @test pdf(PlackettCopula(0.5), [u[i], v[i]]) ≈ l2[i] - end + for i in 1:6 + @test pdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i] + @test pdf(PlackettCopula(0.5), [u[i], v[i]]) ≈ l2[i] end +end - @testset "Sampling" begin - using Random - n_samples = 100 - c = PlackettCopula(0.8) - samples = rand(c, n_samples) - @test size(samples) == (2, n_samples) - end +@testitem "PlackettCopula Sampling" begin + using Random + n_samples = 100 + c = PlackettCopula(0.8) + samples = rand(c, n_samples) + @test size(samples) == (2, n_samples) end From 969ac0ea90a4dfe2f4901b8f45f013aff0b1e8e4 Mon Sep 17 00:00:00 2001 From: Oskar Laverny Date: Thu, 26 Oct 2023 15:11:18 +0200 Subject: [PATCH 5/7] Fix a few typos --- src/Copulas.jl | 2 +- src/MiscellaneousCopulas/PlackettCopula.jl | 14 +++++++------- test/biv_Plackett.jl | 12 +++++++----- 3 files changed, 15 insertions(+), 13 deletions(-) diff --git a/src/Copulas.jl b/src/Copulas.jl index fa2f6e79..b8b18fa0 100644 --- a/src/Copulas.jl +++ b/src/Copulas.jl @@ -62,6 +62,6 @@ end Copulas include("MiscellaneousCopulas/PlackettCopula.jl") export MCopula, WCopula, - SurvivalCopula + SurvivalCopula, PlackettCopula end diff --git a/src/MiscellaneousCopulas/PlackettCopula.jl b/src/MiscellaneousCopulas/PlackettCopula.jl index f174436c..610cfde6 100644 --- a/src/MiscellaneousCopulas/PlackettCopula.jl +++ b/src/MiscellaneousCopulas/PlackettCopula.jl @@ -6,7 +6,7 @@ struct PlackettCopula{P} <: Copula{2} # since it is only bivariate. θ::P # Copula parameter - function PlackettCopula(θ) where {P} + function PlackettCopula(θ) if θ == 1 return IndependentCopula(2) elseif θ == 0 @@ -54,9 +54,9 @@ function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, x::Abstract u = rand(rng) t = rand(rng) a = t * (1 - t) - b = c.θ + a * (c.θ - 1)^2 - cc = 2a * (u * c.θ^2 + 1 - u) + c.θ * (1 - 2a) - d = sqrt(c.θ) * sqrt(c.θ + 4a * u * (1 - u) * (1 - c.θ)^2) + b = C.θ + a * (C.θ - 1)^2 + cc = 2a * (u * C.θ^2 + 1 - u) + C.θ * (1 - 2a) + d = sqrt(C.θ) * sqrt(C.θ + 4a * u * (1 - u) * (1 - C.θ)^2) v = (cc - (1 - 2t) * d) / (2b) x[1] = u x[2] = v @@ -67,9 +67,9 @@ function Base.rand(rng::Distributions.AbstractRNG,C::CT) where CT<: PlackettCopu u = rand(rng) t = rand(rng) a = t * (1 - t) - b = c.θ + a * (c.θ - 1)^2 - cc = 2a * (u * c.θ^2 + 1 - u) + c.θ * (1 - 2a) - d = sqrt(c.θ) * sqrt(c.θ + 4a * u * (1 - u) * (1 - c.θ)^2) + b = C.θ + a * (C.θ - 1)^2 + cc = 2a * (u * C.θ^2 + 1 - u) + C.θ * (1 - 2a) + d = sqrt(C.θ) * sqrt(C.θ + 4a * u * (1 - u) * (1 - C.θ)^2) v = (cc - (1 - 2t) * d) / (2b) x[1] = u x[2] = v diff --git a/test/biv_Plackett.jl b/test/biv_Plackett.jl index c03dc960..ab957f07 100644 --- a/test/biv_Plackett.jl +++ b/test/biv_Plackett.jl @@ -1,11 +1,12 @@ @testitem "PlackettCopula Constructor" begin - @test_broken isa(PlackettCopula(1), Independence) - @test_broken isa(PlackettCopula(Inf),WCopula) # should work in any dimenisons if theta is smaller than the bound. - @test_broken isa(PlackettCopula(0),MCopula) + @test isa(PlackettCopula(1), IndependentCopula) + @test isa(PlackettCopula(Inf),WCopula) # should work in any dimenisons if theta is smaller than the bound. + @test isa(PlackettCopula(0),MCopula) @test_throws ArgumentError PlackettCopula(-0.5) end @testitem "PlackettCopula CDF" begin + using Distributions u = 0.1:0.18:1 v = 0.4:0.1:0.9 l1 = [0.0553778 ,0.1743884 ,0.3166277 ,0.4823228 ,0.6743114 ,0.9 ] @@ -18,6 +19,7 @@ end end @testitem "PlackettCopula PDF" begin + using Distributions u = 0.1:0.18:1 v = 0.4:0.1:0.9 l1 = [1.059211 ,1.023291 ,1.038467 ,1.110077 ,1.272959 ,1.652893 ] @@ -32,7 +34,7 @@ end @testitem "PlackettCopula Sampling" begin using Random n_samples = 100 - c = PlackettCopula(0.8) - samples = rand(c, n_samples) + C = PlackettCopula(0.8) + samples = rand(C, n_samples) @test size(samples) == (2, n_samples) end From f3b3fc2ac56f6ecd6562fd4abd50faad64d7ccb4 Mon Sep 17 00:00:00 2001 From: Oskar Laverny Date: Thu, 26 Oct 2023 16:57:22 +0200 Subject: [PATCH 6/7] make values more precise to pass tests --- test/biv_Plackett.jl | 36 ++++++++++++++++++++++++++++++++---- 1 file changed, 32 insertions(+), 4 deletions(-) diff --git a/test/biv_Plackett.jl b/test/biv_Plackett.jl index ab957f07..a6f6657e 100644 --- a/test/biv_Plackett.jl +++ b/test/biv_Plackett.jl @@ -9,8 +9,22 @@ end using Distributions u = 0.1:0.18:1 v = 0.4:0.1:0.9 - l1 = [0.0553778 ,0.1743884 ,0.3166277 ,0.4823228 ,0.6743114 ,0.9 ] - l2 = [0.02620873, 0.1056116, 0.2349113, 0.4162573, 0.6419255 ,0.9] + l1 = [ + 0.055377800527509735, + 0.1743883734874062, + 0.3166277269195278, + 0.48232275012183223, + 0.6743113969874872, + 0.8999999999999999 + ] + l2 = [ + 0.026208734813001233, + 0.10561162651259381, + 0.23491134194308438, + 0.4162573282722253, + 0.6419254774317229, + 0.9 + ] for i in 1:6 @test cdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i] @@ -22,8 +36,22 @@ end using Distributions u = 0.1:0.18:1 v = 0.4:0.1:0.9 - l1 = [1.059211 ,1.023291 ,1.038467 ,1.110077 ,1.272959 ,1.652893 ] - l2 = [0.8446203 ,1.023291, 1.064891, 0.9360171, 0.7346612, 0.5540166] + l1 = [ + 1.0592107420343486, + 1.023290881054283, + 1.038466936984394, + 1.1100773231007635, + 1.2729591789643138, + 1.652892561983471 + ] + l2 = [ + 0.8446203068160272, + 1.023290881054283, + 1.0648914416282562, + 0.9360170818943749, + 0.7346611825055718, + 0.5540166204986149 + ] for i in 1:6 @test pdf(PlackettCopula(2.0), [u[i], v[i]]) ≈ l1[i] From 9494de852a075d2f5dfb7b6a71fd92a7f4d52bf5 Mon Sep 17 00:00:00 2001 From: Oskar Laverny Date: Thu, 26 Oct 2023 16:57:48 +0200 Subject: [PATCH 7/7] Remove Base.rand --- src/MiscellaneousCopulas/PlackettCopula.jl | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/src/MiscellaneousCopulas/PlackettCopula.jl b/src/MiscellaneousCopulas/PlackettCopula.jl index 610cfde6..02d2b2ec 100644 --- a/src/MiscellaneousCopulas/PlackettCopula.jl +++ b/src/MiscellaneousCopulas/PlackettCopula.jl @@ -62,19 +62,6 @@ function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, x::Abstract x[2] = v return x end -function Base.rand(rng::Distributions.AbstractRNG,C::CT) where CT<: PlackettCopula - x = rand(rng,length(C)) - u = rand(rng) - t = rand(rng) - a = t * (1 - t) - b = C.θ + a * (C.θ - 1)^2 - cc = 2a * (u * C.θ^2 + 1 - u) + C.θ * (1 - 2a) - d = sqrt(C.θ) * sqrt(C.θ + 4a * u * (1 - u) * (1 - C.θ)^2) - v = (cc - (1 - 2t) * d) / (2b) - x[1] = u - x[2] = v - return x -end # Calculate Spearman's rho based on the PlackettCopula parameters function ρ(c::PlackettCopula{P}) where P