[Feat] Add sampling to generic conditionals

src/Conditioning.jl

@@ -157,6 +157,24 @@ end
 function _cdf(CC::ConditionalCopula{d,D,p,T}, v::AbstractVector{<:Real}) where {d,D,p,T}
     return _partial_cdf(CC.C, CC.is, CC.js, Distributions.quantile.(CC.distortions, v), CC.uⱼₛ) / CC.den
 end
+# Sampling: sequential inverse-CDF using conditional distortions
+function Distributions._rand!(rng::Distributions.AbstractRNG, CC::ConditionalCopula{d, D, p}, x::AbstractVector{T}) where {T<:Real, d, D, p}
+    # We want a sample from the COPULA of the conditional model. Let U be a
+    # draw from the conditional joint H_{I|J}(· | u_J). The corresponding
+    # copula coordinates are V_k = F_{i_k|J}(U_k | u_J) = cdf(distortions[k], U_k).
+    # Sample U sequentially by conditioning on J ∪ previously sampled I.
+    J = [j for j in CC.js]
+    ujs =  [u for u in CC.uⱼₛ]
+    for k in 1:d
+        iₖ = CC.is[k]
+        uₖ = rand(rng, DistortionFromCop(CC.C, Tuple(J), Tuple(ujs), iₖ))
+        xₖ = Distributions.cdf(CC.distortions[k], uₖ)
+        x[k] = xₖ
+        push!(J, iₖ)
+        push!(ujs, uₖ)
+    end
+    return x
+end
 
 ###########################################################################
 #####  condition() function

src/Generator.jl

@@ -45,7 +45,12 @@ max_monotony(G::Generator) = throw("This generator does not have a defined max m
 ϕ⁽¹⁾(G::Generator, t) = ForwardDiff.derivative(x -> ϕ(G,x), t)
 ϕ⁻¹⁽¹⁾(G::Generator, t) = ForwardDiff.derivative(x -> ϕ⁻¹(G, x), t)
 ϕ⁽ᵏ⁾(G::Generator, k::Int, t) = taylor(ϕ(G), t, k)[end] * factorial(k)
-ϕ⁽ᵏ⁾⁻¹(G::Generator, k::Int, t; start_at=t) = Roots.find_zero(x -> ϕ⁽ᵏ⁾(G, k, x) - t, start_at)
+ϕ⁽ᵏ⁾⁻¹(G::Generator, k::Int, t; start_at=t) = try 
+    Roots.find_zero(x -> ϕ⁽ᵏ⁾(G, k, x) - t, start_at)
+catch
+    Roots.find_zero(x -> ϕ⁽ᵏ⁾(G, k, x) - t, (0,Inf))
+end
+
 
 
 

test/GenericTests.jl

@@ -610,6 +610,9 @@ Bestiary = filter(GenericTestFilter, Bestiary)
                 i = 3-j
                 for v in (0.3, 0.7)
                     Dd = Copulas.condition(C, j, v)
+                    if !(C isa EmpiricalCopula)
+                        @test all(0 .≤ rand(rng, Dd, 2) .≤ 1) # to ensure the conditional distribution can be sampled. 
+                    end
                     vals = cdf.(Ref(Dd), us)
                     @test all(0.0 .<= vals .<= 1.0)
                     @test all(diff(collect(vals)) .>= -1e-10)
@@ -631,6 +634,9 @@ Bestiary = filter(GenericTestFilter, Bestiary)
             js = tuple(collect(3:d)...)
             ujs = tuple(collect(0.25 + 0.5*rand(rng) for _ in js)...)  # interior values
             CC = condition(C, js, ujs)
+            if !(C isa EmpiricalCopula)
+                @test all(0 .≤ rand(rng, CC, 2) .≤ 1) # to ensure the conditional distribution can be sampled. 
+            end
             pts = [[0.2,0.3], [0.5,0.5], [0.8,0.6]]
             vals = cdf.(CC.C, pts) # only the conditional copula. 
             m_fast = which(Copulas.ConditionalCopula, (CT,                NTuple{d-2, Int}, NTuple{d-2, Float64}))