Add pdf() and logpdf() for Distortions

src/Conditioning.jl

@@ -35,7 +35,7 @@ _process_tuples(::Val{D}, j::Int64, uj::Real) where {D} = ((j,), (uj,))
 function _process_tuples(::Val{D}, js, ujs) where D
     p, p2 = length(js), length(ujs)
     @assert 0 < p < D "js=$(js) must be a non-empty proper subset of 1:D of length at most D-1 (D = $D)"
-    @assert p == p2 && all(0 .<= ujs .<= 1) "uⱼₛ must be in [0,1] and match js length"
+    @assert p == p2 "uⱼₛ length must match js length"
     jst = Tuple(collect(Int, js))
     @assert all(in(1:D), jst)
     ujst = Tuple(collect(float.(ujs)))
@@ -72,6 +72,14 @@ end
 Distributions.logcdf(d::Distortion, t::Real) = log(Distributions.cdf(d, t))
 Distributions.cdf(d::Distortion, t::Real) = exp(Distributions.logcdf(d, t))
 
+# These slow versions are given, but you should probably overrid them: 
+function Distributions.logpdf(d::Distortion, u::Real)
+    (0.0 <= u <= 1.0) || return -Inf
+    v = ForwardDiff.derivative(t -> Distributions.cdf(d, t), float(u))
+    v <= 0 && return -Inf
+    return log(v)
+end
+
 """
     DistortionFromCop{TC,p,T} <: Distortion
 
@@ -127,6 +135,10 @@ struct DistortedDist{Disto, Distrib}<:Distributions.ContinuousUnivariateDistribu
 end
 Distributions.cdf(D::DistortedDist, t::Real) = Distributions.cdf(D.D, Distributions.cdf(D.X, t))
 Distributions.quantile(D::DistortedDist, α::Real) = Distributions.quantile(D.X, Distributions.quantile(D.D, α))
+function Distributions.logpdf(D::DistortedDist, t::Real)
+    u = Distributions.cdf(D.X, t)
+    return Distributions.logpdf(D.X, t) + Distributions.logpdf(D.D, u)
+end
 
 """
     ConditionalCopula{d} <: Copula{d}

src/SklarDist.jl

@@ -58,6 +58,7 @@ SklarDist(C, m) = SklarDist(C, Tuple(m))
 Base.length(S::SklarDist{CT,TplMargins}) where {CT,TplMargins} = length(S.C)
 Base.eltype(S::SklarDist{CT,TplMargins}) where {CT,TplMargins} = Base.eltype(S.C)
 Distributions.cdf(S::SklarDist{CT,TplMargins},x) where {CT,TplMargins} = Distributions.cdf(S.C, [Distributions.cdf(mi, xi) for (mi, xi) in zip(S.m, x)])
+Distributions.logcdf(S::SklarDist{CT,TplMargins},x) where {CT,TplMargins} = log(Distributions.cdf(S, x))
 function Distributions._rand!(rng::Distributions.AbstractRNG, S::SklarDist{CT,TplMargins}, x::AbstractVector{T}) where {CT,TplMargins,T}
     Random.rand!(rng,S.C,x)
     clamp!(x, nextfloat(T(0)), prevfloat(T(1)))

src/UnivariateDistribution/Distortions/ArchimedeanDistortion.jl

@@ -16,3 +16,8 @@ function Distributions.quantile(D::ArchimedeanDistortion{TG, T}, α::Real) where
     return ϕ(D.G, y - D.sJ)
 end
 ## ConditionalCopula moved next to ArchimedeanCopula definition
+function Distributions.logpdf(D::ArchimedeanDistortion{TG, T}, u::Real) where {TG, T}
+    ξ = ϕ⁻¹(D.G, float(u))
+    num = ϕ⁽ᵏ⁾(D.G, D.p + 1, D.sJ + ξ)
+    return log(abs(num)) - log(abs(D.den)) - log(abs(ϕ⁽ᵏ⁾(D.G, 1, ξ)))
+end

src/UnivariateDistribution/Distortions/BivArchimaxDistortion.jl

@@ -36,3 +36,38 @@ function Distributions.logcdf(D::BivArchimaxDistortion, z::Real)
     r = max(r, T(0))
     return min(log(-ϕ⁽¹⁾(D.gen, S * A0)) + log(-ϕ⁻¹⁽¹⁾(D.gen, D.uⱼ) * r), T(0))
 end
+function Distributions.logpdf(D::BivArchimaxDistortion, z::Real)
+    T = typeof(z)
+    z ≤ 0    && return T(-Inf)
+    z ≥ 1    && return T(-Inf)
+    D.uⱼ ≤ 0 && return T(-Inf)
+    D.uⱼ ≥ 1 && return T(0)
+
+    x = ϕ⁻¹(D.gen, z)
+    y = ϕ⁻¹(D.gen, D.uⱼ)
+    S = x + y
+    S <= 0 && return T(-Inf)
+    t  = D.j==2 ? _safett(y / S) : _safett(x / S)
+
+    A0 = A(D.tail, t)
+    A1 = dA(D.tail, t)
+    A2 = d²A(D.tail, t)
+    r = D.j==2 ? (A0 + (1 - t) * A1)  : (A0 - t * A1)
+    r = max(r, T(0))
+
+    G = S * A0
+    phi1G = ϕ⁽¹⁾(D.gen, G)
+    phi2G = ϕ⁽ᵏ⁾(D.gen, 2, G)
+
+    inv_cond = ϕ⁻¹⁽¹⁾(D.gen, D.uⱼ)   # derivative of inverse at conditioning value
+    dx_dz   = ϕ⁻¹⁽¹⁾(D.gen, z)      # derivative of inverse at z
+
+    dGdx = D.j==2 ? (A0 - t * A1) : (A0 + (1 - t) * A1)
+    drdx = - (t * (1 - t) / S) * A2
+
+    inner = phi2G * dGdx * r + phi1G * drdx
+    val = inv_cond * dx_dz * inner
+
+    val <= zero(T) && return T(-Inf)
+    return T(log(val))
+end

src/UnivariateDistribution/Distortions/BivEVDistortion.jl

@@ -34,4 +34,54 @@ function Distributions.logcdf(D::BivEVDistortion, z::Real)
 
     # upper clip but no lower clip
     return min(logval + log(max(tolog, T(0))), T(0))
+end
+function Distributions.logpdf(D::BivEVDistortion, z::Real)
+    T = typeof(z)
+    # bounds and degeneracies
+    z ≤ 0    && return T(-Inf)
+    z ≥ 1    && return T(-Inf)
+    D.uⱼ ≤ 0 && return T(-Inf)
+    D.uⱼ ≥ 1 && return T(0)  # conditioning on 1 ⇒ uniform density = 1 => logpdf = 0
+
+    if D.j == 2
+        # Condition on the second variable : V = D.uⱼ, free = u=z
+        x, y = -log(z), -log(D.uⱼ)
+        s = x + y
+        w = x / s
+        Aw, dAw = A(D.tail, w), dA(D.tail, w)
+        ddAw = d²A(D.tail, w)
+        Tval = Aw - w * dAw
+        Tval ≤ 0 && return T(-Inf)
+
+        logval = -s * Aw + y
+        # derivatives
+        # logval' = (Aw + (y/s)*dAw) / z
+        lvp = (Aw + (y / s) * dAw) / z
+        # T'(w)*dw/dz = w * ddAw * y / (z * s^2)
+        tp_term = w * ddAw * y / (z * s^2)
+        B = tp_term + Tval * lvp
+        B ≤ 0 && return T(-Inf)
+
+        return logval + log(B)
+    else
+        # Condition on the first variable : U = D.uⱼ, free = v=z
+        x, y = -log(D.uⱼ), -log(z)
+        s = x + y
+        w = x / s
+        Aw, dAw = A(D.tail, w), dA(D.tail, w)
+        ddAw = d²A(D.tail, w)
+        Tval = Aw + (1 - w) * dAw
+        Tval ≤ 0 && return T(-Inf)
+
+        logval = -s * Aw + x
+        # derivatives
+        # logval' = (Aw - (x/s)*dAw) / z
+        lvp = (Aw - (x / s) * dAw) / z
+        # T'(w)*dw/dz = x * (1 - w) * ddAw / (z * s^2)
+        tp_term = x * (1 - w) * ddAw / (z * s^2)
+        B = tp_term + Tval * lvp
+        B ≤ 0 && return T(-Inf)
+
+        return logval + log(B)
+    end
 end

src/UnivariateDistribution/Distortions/BivFGMDistortion.jl

@@ -16,3 +16,9 @@ function Distributions.quantile(D::BivFGMDistortion, α::Real)
     # Solve a u^2 - (1+a) u + α = 0 and pick the root in [0,1]
     return ((1 + a) - sqrt((1 + a)^2 - 4*a*α)) / 2a
 end
+function Distributions.logpdf(D::BivFGMDistortion, u::Real)
+    v = D.θ * (1 - 2 * D.uⱼ) * (1 - 2 * u)
+    p = 1 + v
+    p <= 0 && return -Inf
+    return log1p(v)
+end

src/UnivariateDistribution/Distortions/FlipDistortion.jl

@@ -8,3 +8,4 @@ Distributions.cdf(D::FlipDistortion, u::Real) = 1.0 - Distributions.cdf(D.base,
 Distributions.quantile(D::FlipDistortion, α::Real) = 1.0 - Distributions.quantile(D.base, 1.0 - float(α))
 
 ## Methods moved next to SurvivalCopula type
+Distributions.logpdf(D::FlipDistortion, u::Real) = Distributions.logpdf(D.base, 1.0 - float(u))

src/UnivariateDistribution/Distortions/GaussianDistortion.jl

@@ -20,3 +20,8 @@ function (D::GaussianDistortion)(X::Distributions.Normal)
     return Distributions.Normal(μ + σ*D.μz, σ*D.σz)
 end
 ## Methods moved to EllipticalCopulas/GaussianCopula.jl
+function Distributions.logpdf(d::GaussianDistortion, u::Real)
+    N = Distributions.Normal()
+    q = Distributions.quantile(N, u)
+    return Distributions.logpdf(N, (q - d.μz)/d.σz) - log(abs(d.σz)) - Distributions.logpdf(N, q)
+end

src/UnivariateDistribution/Distortions/HistogramDistortion.jl

@@ -44,3 +44,18 @@ end
     k = idx - 1
     return (k + frac) / d.m
 end
+@inline function Distributions.logpdf(d::HistogramBinDistortion, u::Real)
+    u_f = float(u)
+    if !isfinite(u_f) || u_f < 0.0 || u_f > 1.0
+        return -Inf
+    end
+    m = d.m
+    s = m * u_f
+    k = min(max(floor(Int, s), 0), m - 1)
+    idx = k + 1
+    p = d.probs[idx]
+    if p <= 0.0
+        return -Inf
+    end
+    return log(m) + log(float(p))
+end

src/UnivariateDistribution/Distortions/MDistortion.jl

@@ -7,3 +7,4 @@ struct MDistortion{T} <: Distortion
 end
 Distributions.cdf(D::MDistortion, u::Real) = min(u / D.v, 1)
 Distributions.quantile(D::MDistortion, α::Real) = α * D.v
+Distributions.logpdf(D::MDistortion, u::Real) = (0 <= u <= D.v) ? -log(D.v) : -Inf

src/UnivariateDistribution/Distortions/NoDistortion.jl

@@ -6,3 +6,4 @@ Distributions.cdf(::NoDistortion, u::Real) = u
 Distributions.quantile(::NoDistortion, α::Real) = α
 (::NoDistortion)(X::Distributions.UnivariateDistribution) = X
 (::NoDistortion)(::Distributions.Uniform) = Distributions.Uniform()
+Distributions.logpdf(::NoDistortion, u::Real) = zero(u)

src/UnivariateDistribution/Distortions/PlackettDistortion.jl

@@ -21,3 +21,29 @@ function Distributions.logcdf(D::PlackettDistortion, u::Real)
     return num - den
 end
 ## DistortionFromCop moved next to PlackettCopula
+function Distributions.logpdf(D::PlackettDistortion, u::Real)
+    θ, v = D.θ, D.uⱼ
+    η = θ - one(θ)
+
+    # θ == 1 (independence) -> uniform conditional density
+    if η == zero(η)
+        return zero(eltype(θ))
+    end
+
+    t1 = η * (u + v) + one(θ)
+    s1 = sqrt(max(zero(θ), t1 * t1 - 4θ * η * u * v))
+
+    # A(u) = 1 - (t1 - 2θ*u)/s1
+    # derivative A'(u) = -d/du[(t1 - 2θ*u)/s1]
+    dt1 = η
+    ds1 = η * (t1 - 2θ * v) / s1
+    dB = ((dt1 - 2θ) * s1 - (t1 - 2θ * u) * ds1) / (s1 * s1)
+    Ap = -dB
+
+    t2 = η * (1 + v) + one(θ)
+    s2 = sqrt(max(zero(θ), t2 * t2 - 4θ * η * v))
+
+    denomA = 1 - (t2 - 2θ) / s2
+
+    return log(abs(Ap)) - log(denomA)
+end

src/UnivariateDistribution/Distortions/StudentDistortion.jl

@@ -18,3 +18,9 @@ function Distributions.quantile(d::StudentDistortion, α::Real)
     return Distributions.cdf(Tu, d.μz + d.σz * zα)
 end
 ## Methods moved next to TCopula type
+function Distributions.logpdf(d::StudentDistortion, u::Real)
+    Tu = Distributions.TDist(d.ν)
+    Tcond = Distributions.TDist(d.νp)
+    z = Distributions.quantile(Tu, float(u))
+    return Distributions.logpdf(Tcond, (z - d.μz) / d.σz) - log(abs(d.σz)) - Distributions.logpdf(Tu, z)
+end

src/UnivariateDistribution/Distortions/WDistortion.jl

@@ -7,3 +7,4 @@ struct WDistortion{T} <: Distortion
 end
 Distributions.cdf(D::WDistortion, u::Real) = max(u + D.v - 1, 0) / D.v
 Distributions.quantile(D::WDistortion, α::Real) = α * D.v + (1 - D.v)
+Distributions.logpdf(D::WDistortion, u::Real) = (u + D.v > 1 ? -log(D.v) : -Inf)

test/GenericTests.jl

@@ -614,13 +614,20 @@ Bestiary = filter(GenericTestFilter, Bestiary)
                         @test all(0 .≤ rand(rng, Dd, 2) .≤ 1) # to ensure the conditional distribution can be sampled. 
                     end
                     vals = cdf.(Ref(Dd), us)
+                    pvals = pdf.(Ref(Dd), us)
                     @test all(0.0 .<= vals .<= 1.0)
                     @test all(diff(collect(vals)) .>= -1e-10)
+                    @test all(pvals .>= 0)
                     if check_biv_conditioning(C) && has_spec
                         Dgen  = @invoke Copulas.DistortionFromCop(C::Copulas.Copula{d}, (j,), (v,), i)
                         vals_gen  = cdf.(Ref(Dgen),  us)
+                        pvals_gen = pdf.(Ref(Dgen), us)
+                        tol = C isa Copulas.GaussianCopula ? 1e-2 : 1e-3
                         for (vf, vg) in zip(vals, vals_gen)
-                            @test isapprox(vf, vg, atol=1e-3, rtol=1e-3)
+                            @test isapprox(vf, vg, atol=tol, rtol=tol)
+                        end
+                        for (vf, vg) in zip(pvals, pvals_gen)
+                            @test isapprox(vf, vg, atol=tol, rtol=tol)
                         end
                     elseif CT <: Copulas.MCopula
                         @test all(vals .≈ min.(collect(us) ./ v, 1))