diff --git a/src/ArchimedeanCopula.jl b/src/ArchimedeanCopula.jl index e0515a61..51234dbf 100644 --- a/src/ArchimedeanCopula.jl +++ b/src/ArchimedeanCopula.jl @@ -66,26 +66,18 @@ function Distributions._logpdf(C::ArchimedeanCopula{d,TG}, u) where {d,TG} if !all(0 .<= u .<= 1) return eltype(u)(-Inf) end - sum_ϕ⁻¹u = 0.0 - sum_logϕ⁽¹⁾ϕ⁻¹u = 0.0 + T = promote_type(Float64, eltype(u)) # the FLoat64 here should be eltype(C) when copulas wil be type agnostic... + logdenom = sum_ϕ⁻¹u = zero(T) for us in u ϕ⁻¹u = ϕ⁻¹(C,us) sum_ϕ⁻¹u += ϕ⁻¹u - sum_logϕ⁽¹⁾ϕ⁻¹u += log(-ϕ⁽¹⁾(C,ϕ⁻¹u)) # log of negative here because ϕ⁽¹⁾ is necessarily negative + logdenom += log(-ϕ⁽¹⁾(C,ϕ⁻¹u)) # log of negative here because ϕ⁽¹⁾ is necessarily negative end - numer = ϕ⁽ᵏ⁾(C, d, sum_ϕ⁻¹u) - dimension_sign = iseven(d) ? 1.0 : -1.0 #need this for log since (-1.0)ᵈ ϕ⁽ᵈ⁾ ≥ 0.0 - - - # I am not sure this is the right reasoning : - if numer == 0 - if sum_logϕ⁽¹⁾ϕ⁻¹u == -Inf - return Inf - else - return -Inf - end + numer = abs(ϕ⁽ᵏ⁾(C, d, sum_ϕ⁻¹u)) + if numer > 0 + return log(numer) - logdenom else - return log(dimension_sign*numer) - sum_logϕ⁽¹⁾ϕ⁻¹u + return -T(Inf) end end diff --git a/test/ClaytonCopula.jl b/test/ClaytonCopula.jl index e36c8dc6..a4c39862 100644 --- a/test/ClaytonCopula.jl +++ b/test/ClaytonCopula.jl @@ -6,8 +6,8 @@ y = x cdf1 = [0.0, 0.1796053020267749, 0.37796447300922725, 0.6255432421712244, 1.0] cdf2 = [0.0, 0.0, 0.17157287525381, 0.5358983848622453, 1.0] - pdf1 = [Inf, 2.2965556205046926, 1.481003649342278, 1.614508582188617, 3.0] - pdf2 = [Inf, 0.0, 1.0, 2 / 3, 0.5] + pdf1 = [0.0, 2.2965556205046926, 1.481003649342278, 1.614508582188617, 3.0] + pdf2 = [0.0, 0.0, 1.0, 2 / 3, 0.5] for i in 1:5 @test cdf(ClaytonCopula(2,2),[x[i],y[i]]) ≈ cdf1[i] @test cdf(ClaytonCopula(2,-0.5),[x[i],y[i]]) ≈ cdf2[i]