# Issues implementing Copula regression using copula and Turing

**URL:** <https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907>\
**Category:** Probabilistic Programming\
**Tags:** question, turing, copula\
**Created:** [November 27, 2022, 3:45pm UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907 "2022-11-27T15:45:14Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Kokora](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kokora/32/7867_2.png) [@Kokora](https://discourse.julialang.org/u/Kokora)\
**Post date:** [November 27, 2022, 3:45pm UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/1 "2022-11-27T15:45:15Z")

</div>

I am trying to implement a copula regression model by leveraging Turing and copulas in Julia.

The author of copula implemented a nice working example as below

```julia
using Copulas, Distributions, Random, Turing, Plots, StatsPlots

Random.seed!(123)
D = SklarDist(ClaytonCopula(3,7), (Exponential(1.0), Pareto(3.0), Exponential(2.0)))
draws = rand(D, 2_000)

@model function copula(X)
    # Priors
    θ ~ TruncatedNormal(1.0, 1.0, 0, Inf)
    θ₁ ~ TruncatedNormal(1.0, 1.0, 0, Inf)
    θ₂ ~ TruncatedNormal(1.0, 1.0, 0, Inf)
    θ₃ ~ TruncatedNormal(1.0, 1.0, 0, Inf)

    # Marginals distributions and copula
    X₁ = Exponential(θ₁)
    X₂ = Pareto(θ₂)
    X₃ = Exponential(θ₃)
    C = ClaytonCopula(3,θ)
    D = SklarDist(C, (X₁, X₂, X₃))
    Turing.Turing.@addlogprob! loglikelihood(D, X)
end

sampler = NUTS() # MH() works too
chain = sample(copula(draws), sampler, MCMCThreads(), 1_00, 4)
plot(chain)

```

I would like to implement a regression equivalent of that model. I have a simple example below. But NUTS and HMC samplers do not work. PG works but does not produce a reliable results.

```julia
using Random
Random.seed!(43);

n = 200
p = 2

X = randn(n, p)
β1 = randn(p)

Y = X*hcat(β1,β1) + transpose(rand(SklarDist(GaussianCopula([1 .8; .8 1]), (Normal(0.0, 1.), Normal(0.0,1.0))), n))

#-- Look at the data
cornerplot(Y, compact=true)

#-- Model
@model function copulaReg(Y, X)

    n, K = size(X)

    # Priors
    ρ ~ truncated(Normal(0,.3), -1, 1) #-- Corelated
    θ ~ Normal(0.,1.)
    θ₁ ~ Normal(0.,1.)
    θ₂ ~ Normal(0.,1.)

    β1 ~ filldist(Normal(0., 1.), K)
    β2 ~ filldist(Normal(0., 1.), K)
     β ~ filldist(Normal(0., 1.), K)

    #-- Regression parameters
    μ1 = zeros(n)
    μ2 = zeros(n)

    for i in 1:n
        μ1[i] = θ₁ .+ sum(X[i,:].* β1)
        μ2[i] = θ₂ .+ sum(X[i,:].* β2)

        #Turing.Turing.@addlogprob! loglikelihood(SklarDist(GaussianCopula([1. ρ; ρ 1.0]), (Normal(μ1[i], σ1), Normal(μ2[i], σ2))), Y[i,:])
        Turing.Turing.@addlogprob! loglikelihood(SklarDist(GaussianCopula([1. ρ; ρ 1.0]), (Normal(μ1[i], 1.), Normal(μ2[i], 1.))), Y[i,:])
    end
end

#-- sampler
sampler = NUTS() # doe not work
sampler = PG(10) # works but gives weird results!!
#-- sample from the posterior
chain = sample(copulaReg(Y, X), sampler, MCMCThreads(), 2_000, 3)

describe(chain)

```

I am not sure why NUTS and HMC are not working on this. Any help on this will be greatly appreciated. Thank you.

---

<div class="post-metadata">

**Author:** ![bertschi](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bertschi/32/33462_2.png) [@bertschi](https://discourse.julialang.org/u/bertschi)\
**Post date:** [November 27, 2022, 7:22pm UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/2 "2022-11-27T19:22:08Z")

</div>

Your code is not generic enough, i.e., `zeros` creates an array of `Float64` and thus fails when Turing tries to take gradients using ForwardDiff by default as `ForwardDiff.Dual` cannot be written to those float arrays. Try with `μ1 = zeros(typeof(θ₁), n)` (same for `μ2`) or don’t use these arrays at all, but simply use local variables within the loop, i.e., `μ1 = θ₁ .+ sum(X[i,:].* β1)` (same for `μ2`).  
I also needed to fix the correlation matrix of the Gaussian copula to `GaussianCopula([1+eps ρ; ρ 1+eps])` for some small `eps`, e.g., `1e-7`, in order to fix numerical issues when `ρ` gets close to one.

---

<div class="post-metadata">

**Author:** ![Kokora](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kokora/32/7867_2.png) [@Kokora](https://discourse.julialang.org/u/Kokora)\
**Post date:** [November 27, 2022, 10:43pm UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/3 "2022-11-27T22:43:45Z")

</div>

> [@bertschi](#):
>
> GaussianCopula([1+eps ρ; ρ 1+eps])

@bertschi haha - I see now! That works just perfect! Thank you.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [November 28, 2022, 8:57am UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/4 "2022-11-28T08:57:48Z")

</div>

@kokora, if I may what was the purpose of this MWE ?

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [November 28, 2022, 10:02am UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/5 "2022-11-28T10:02:09Z")

</div>

The trick with `eps` works great. Do you think changing `ρ` to `ρ ~ truncated(Normal(0,.3), -0.99, 0.99)` would also work? In some sense this is not an ‘uninformed’ prior, but if the posterior bangs against the truncation, it would indicate the model needs adjustment.  
I think having diagonal \> 1 on correlation matrix feels too tricky.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [November 28, 2022, 10:07am UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/6 "2022-11-28T10:07:42Z")

</div>

> [@Dan](#):
>
> I think having diagonal \> 1 on correlation matrix feels too tricky.

The matrix you input is re-normalized to be a propper correlation matrix by the copula constructor. Thus, this `+eps` really is a contraction of `rho` by some amount, and it is not great as the variable rho does not mean anything anymore.

Truncating `rho`’s prior, if it works, would be indeed more sensible as the value obtained for it would still be interpretable.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [November 28, 2022, 10:15am UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/7 "2022-11-28T10:15:25Z")

</div>

> [@lrnv](#):
>
> The matrix you input is re-normalized to be a propper correlation matrix by the copula constructor

Perhaps an error/exception is better. This thread might be the example to show why an error/exception for a non-proper correlation matrix is better. The ‘projection’ to a correlation matrix can be done outside the copula constructor (maybe even a CorrelationMatrix type is needed).

Relevant link: [https://en.wikipedia.org/wiki/Correlation#Nearest\_valid\_correlation\_matrix](https://en.wikipedia.org/wiki/Correlation#Nearest_valid_correlation_matrix)

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [November 28, 2022, 10:20am UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/8 "2022-11-28T10:20:38Z")

</div>

You might be right.

Moreover, with `NUTS()`, even if resulsts look great, I still got a bunch of warnings :

```julia
┌ Warning: The current proposal will be rejected due to numerical error(s).
│ isfinite.((θ, r, ℓπ, ℓκ)) = (true, false, false, false)      
└ @ AdvancedHMC C:\Users\lrnv\.julia\packages\AdvancedHMC\iWHPQ\src\hamiltonian.jl:47

```

Which do not look great. What should i do about these ?

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [November 28, 2022, 10:23am UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/9 "2022-11-28T10:23:19Z")

</div>

Got those too. Didn’t dig enough to understand either.

---

<div class="post-metadata">

**Author:** ![Kokora](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kokora/32/7867_2.png) [@Kokora](https://discourse.julialang.org/u/Kokora)\
**Post date:** [November 28, 2022, 10:39pm UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/10 "2022-11-28T22:39:36Z")

</div>

@Dan Good point! I can also constrain that prior away from -1 and 1.  
I am also thinking about specifying a regression model for ρ so that ρ = tanh(XB)\*0.99 with prior on B instead. That 0.99 allows us to stay away from -1 and 1.

---

<div class="post-metadata">

**Author:** ![Kokora](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kokora/32/7867_2.png) [@Kokora](https://discourse.julialang.org/u/Kokora)\
**Post date:** [November 28, 2022, 10:42pm UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/11 "2022-11-28T22:42:25Z")

</div>

@lrnv I conquer - It is not at all clear what these warnings message are about? Maybe I should reopen this and mark as not fully resolved yet?

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [November 28, 2022, 10:47pm UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/12 "2022-11-28T22:47:18Z")

</div>

@Kokora About the warnings I cannot help you, but maybe someone else can. On the other hand, about the prior on `ρ`, you might also try other copula families than the Gaussian on your real data.

---

<div class="post-metadata">

**Author:** ![Dan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dan/32/42581_2.png) [@Dan](https://discourse.julialang.org/u/Dan)\
**Post date:** [November 28, 2022, 11:26pm UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/13 "2022-11-28T23:26:05Z")

</div>

> [@Kokora](#):
>
> ρ = tanh(XB)\*0.99

If using `tanh` formula, it is not necessary perhaps to limit to `0.99` as the problematic value is only reached ‘at infinity’. It also seems `tanh` is being used as a bijector from (-Inf,Inf) to (-1,1). So it might be interesting to study Bijector package to formally make this into one. It might somehow help some gradient calculations.  
[https://github.com/TuringLang/Bijectors.jl](https://github.com/TuringLang/Bijectors.jl)

(I’m learning as I go too).

UPDATE (after a bit of aforementioned ‘learning’):  
`Bijectors.jl` already provides just such a ‘link’ (statistics name for this invertible function). So:

```julia
using Bijectors
cb = Bijectors.CorrBijector()

```

And now:

```julia
julia> Inverse(cb)([0 0.4; 0 0])
2×2 Matrix{Float64}:
 1.0 0.379949
 0.379949 1.0

```

This bijector is good for larger matrices. Basically, the `0.4` in the matrix is the `XB` and the exact `tanh` transform is done. For larger correlation matrices, bigger upper triangular matrices can be supplied.

---

<div class="post-metadata">

**Author:** ![Kokora](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kokora/32/7867_2.png) [@Kokora](https://discourse.julialang.org/u/Kokora)\
**Post date:** [November 29, 2022, 3:55am UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/14 "2022-11-29T03:55:44Z")

</div>

@lrnv Yes, I plan to implement other Copula. I wanted to implement the following models for the rotated Clayton copulas defined as follows:

C90(u1, u2) = u2 − C(1 − u1, u2),  
C180(u1, u2) = u1 + u2 − 1 + C(1 − u1, 1 − u2),  
C270(u1, u2) = u1 − C(u1, 1 − u2),

where C is the standard Clayton Copula.  
See this paper [here](https://doi.org/10.1002/sim.7985)  
These are defined on page 6 of the paper.  
I think these will take some thinking to do!  
Any idea?

---

<div class="post-metadata">

**Author:** ![Kokora](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kokora/32/7867_2.png) [@Kokora](https://discourse.julialang.org/u/Kokora)\
**Post date:** [November 29, 2022, 4:32am UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/15 "2022-11-29T04:32:45Z")

</div>

> [@Dan](#):
>
> ```julia
> using Bijectors
> cb = Bijectors.CorrBijector()
> 
> ```

Ok I see your point. It took a little to get it. But it does not seem to work. See a small example below

```julia
using Copulas, Distributions, Random, Turing, Plots, StatsPlots
n = 200
p = 2

X = randn(n, p)
β1 = randn(p)

Y = X*hcat(β1,β1) + transpose(rand(SklarDist(GaussianCopula([1 .8; .8 1]), (Normal(0.0, 1.), Normal(0.0,1.0))), n))#transpose(rand(MultivariateNormal(zeros(2), [1 0.8; 0.8 1]),n))

#-- Model without bijector for Tanh
@model function copulaReg3(Y, X)

    n, K = size(X)
    cb = Bijectors.CorrBijector()

    # Priors
    #ρ ~ truncated(Normal(0,.3), -.99, 0.99) #-- Corelated
    θ ~ Normal(0.,1.)
    θ₁ ~ Normal(0.,1.)
    θ₂ ~ Normal(0.,1.)

    β1 ~ filldist(Normal(0., 1.), K)
    β2 ~ filldist(Normal(0., 1.), K)
     β ~ filldist(Normal(0., 1.), K)

    #-- Regression parameters
    μ1 = zeros(typeof(θ₁), n)
    μ2 = zeros(typeof(θ₂), n)
    μ3 = zeros(typeof(θ₂), n)
     ρ = zeros(typeof(θ), n)

    for i in 1:n
        μ1[i] = θ₁+sum(X[i,:].* β1)
        μ2[i] = θ₂+sum(X[i,:].* β2)
        μ3[i] = θ+sum(X[i,:].* β)
        ρ[i] = Inverse(cb)([0. μ3[i]; 0. 0.])[2,1]

        #Turing.Turing.@addlogprob! loglikelihood(SklarDist(GaussianCopula([1. ρ; ρ 1.0]), (Normal(μ1[i], σ1), Normal(μ2[i], σ2))), Y[i,:])
        Turing.Turing.@addlogprob! loglikelihood(SklarDist(GaussianCopula([1 ρ[i]; ρ[i] 1]), (Normal(μ1[i], 1.), Normal(μ2[i], 1.))), Y[i,:])
    end
end

#-- sampler
sampler = NUTS()

#-- sample from the posterior
chain = sample(copulaReg3(Y, X), sampler, MCMCThreads(), 2_000, 3)

```

unless there is something I am missing.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [November 29, 2022, 10:34am UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/16 "2022-11-29T10:34:42Z")

</div>

> [@Kokora](#):
>
> Yes, I plan to implement other Copula. I wanted to implement the following models for the rotated Clayton copulas

For that, you might use the `SurvivalCopula` class, that is meant to flip directions of any copula:

```julia
using Copulas
C = ClaytonCopula(2,3.0) # bivariate clayton with theta = 3.0
C90 = SurvivalCopula(C,(1,)) # flips the first dimension
C270 = SurvivalCopula(C,(2,)) # flips only the second dimension. 
C180 = SurvivalCopula(C,(1,2)) # flips both dimensions.

```

You might want to update to Copulas.jl v0.1.7 for that to work correctly.

---

<div class="post-metadata">

**Author:** ![Kokora](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kokora/32/7867_2.png) [@Kokora](https://discourse.julialang.org/u/Kokora)\
**Post date:** [November 29, 2022, 9:31pm UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/17 "2022-11-29T21:31:43Z")

</div>

@lrnv Exactly what I needed. It works just nicely! Thank you.

---

<div class="post-metadata">

**Author:** ![Kokora](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kokora/32/7867_2.png) [@Kokora](https://discourse.julialang.org/u/Kokora)\
**Post date:** [December 9, 2022, 5:26am UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/18 "2022-12-09T05:26:58Z")

</div>

One question - does Copula work with OrderedLogistic distribution?

I have this small example below

```julia
using Copulas, Distributions, Random, Turing, Plots, StatsPlots
using DataFrames
#using CairoMakie
using AlgebraOfGraphics
using Distributions
using StatsFuns#: logit
using Bijectors
using LazyArrays
using LinearAlgebra
using Random

#--- I am trying to fit the following copula model
#-- Parameters 
seed!(123)
n = 500 # sample sizes
p = 3 # Number of caovariates 
β = randn(p)
X = hcat(randn(n, p-1), zeros(n, 1))
lp = X*β
#-- Simulate data from a Clayton copula model

C180 = SurvivalCopula(ClaytonCopula(2,3.0),(1,2)) #ClaytonCopula(2,3.0) # bivariate clayton with theta = 3.0
#D = SklarDist(C180, (Normal(0.0, 1.), Normal(0.0,1.0)))
#D180 = SklarDist(C180, (Normal(0.0, 1.), Normal(0.0,1.0)))

Y = zeros(n, 2) ## Store the data
for i in 1:n
    Y[i,:] = rand( SklarDist(C180, (Normal(lp[i], 1.), Logistic(lp[i],1.0))), 1)
end

#-- Plot the data
cornerplot( Y, compact=true)

#--- Copula model for this joint data
#-- Clayton Copula Model 180 Rotate
@model function Model1Clay180(Y, X)

    n, K = size(X)

    # Priors
    θ ~ Normal(0.,1.)
    θ₁ ~ Normal(0.,1.)
    θ₂ ~ Normal(0.,1.)

    β1 ~ filldist(Normal(0., 1.), K)
    β2 ~ filldist(Normal(0., 1.), K)
     β ~ filldist(Normal(0., 1.), K)

    # Regression parameters
    μ1 = zeros(typeof(θ₁), n)
    μ2 = zeros(typeof(θ₂), n)
     ρ = zeros(typeof(θ), n)

    # Implement Shrinkage Priors
    for i in 1:n
        μ1[i] = θ₁ .+ sum(X[i,:].* β1)
        μ2[i] = θ₂ .+ sum(X[i,:].* β2)
        ρ[i] = exp(θ .+ sum(X[i,:].* β))

        Turing.Turing.@addlogprob! loglikelihood(SklarDist(SurvivalCopula(ClaytonCopula(2,ρ[i]),(1,2)), (Normal(μ1[i], 1.), Logistic(μ2[i], 1.))), Y[i,:])
       
    end
end

# Fit the model to the data
sampler = NUTS()

# sample from the posterior

Model1chain = sample(Model1Clay180(Y, X), sampler, MCMCThreads(), 2_000, 3)

# Second par of the problem
# SUppose Y[:,2] is not full observed but we only observe categories obtained by dichotomizing Y2
using CategoricalArrays

Ycat = cut(Y[:,2], [-Inf, -1.5, 0, 1.5, Inf], labels = [1,2,3,4])

#--- Our model becomes (using OrderedLogistic distribition)

#-- Clayton Copula Model 180 Rotate
@model function Model2Clay180(Y1, Y2, X)

    n, K = size(X)
    ncateg = maximum(Y2)
     
    # Priors
    cutpoints ~ Bijectors.ordered(filldist(Normal(0., 10.), ncateg - 1))
    θ ~ Normal(0.,1.)
    θ₁ ~ Normal(0.,1.)
    θ₂ ~ Normal(0.,1.)

    β1 ~ filldist(Normal(0., 1.), K)
    β2 ~ filldist(Normal(0., 1.), K)
     β ~ filldist(Normal(0., 1.), K)

    # Regression parameters
    μ1 = zeros(typeof(θ₁), n)
    μ2 = zeros(typeof(θ₂), n)
     ρ = zeros(typeof(θ), n)

    # Implement Shrinkage Priors
    for i in 1:n
        μ1[i] = θ₁ .+ sum(X[i,:].* β1)
        μ2[i] = θ₂ .+ sum(X[i,:].* β2)
        ρ[i] = exp(θ .+ sum(X[i,:].* β))

        Turing.Turing.@addlogprob! loglikelihood(SklarDist(SurvivalCopula(ClaytonCopula(2,ρ[i]),(1,2)), (Normal(μ1[i], 1.), OrderedLogistic(μ2[i],cutpoints))), [Y1[i],Y2[i]] )
       
    end
end

#-- Fit the model
Model2chain = sample(Model2Clay180(Y[:,1], Ycat, X), sampler, MCMCThreads(), 2_000, 3)

```

This does not at all for me.

---

<div class="post-metadata">

**Author:** ![lrnv](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/lrnv/32/19373_2.png) [@lrnv](https://discourse.julialang.org/u/lrnv)\
**Post date:** [December 9, 2022, 10:52am UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/19 "2022-12-09T10:52:35Z")

</div>

Hum…

`OrderedLogistic` seem to come from `Turing.jl`, but the cdf corresponding to this distribution is not implemented. We could add it ourselves, but that’ll be type piracy:

```julia
function Distributions.cdf(d::OrderedLogistic{T, Vector{T2}}, x::T3) where {T<:Real, T2<:Real, T3<:Real}
    rez = pdf(d,x)
    for atom in support(d)
        if atom < x
            rez += pdf(d,atom)
        end
    end
    return rez
end

```

But then I still get numerical errors from `Turing.unsafe_logpdf_ordered_logistic`…

I can assure you that this _should_ work, but this distribution does not implement the interface from [Create New Samplers and Distributions · Distributions.jl](https://juliastats.org/Distributions.jl/stable/extends/#Univariate-Distribution) correctly, namely the `cdf` function is not implemented. Maybe we could ping `Turing.jl` with a github issue to fix this.

Another option is to go to `Distributions.jl` and ask for a generic fallback for the cdf of a discrete distribution knowing the pdf, which seems not to be implemented right now (i dont know why).

EDIT: As a side note, your model takes ages to run on my machine. I am not familiar enough with Turing to know if it is normal, but maybe there are things to do to make it a lot faster ?

---

<div class="post-metadata">

**Author:** ![Kokora](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/kokora/32/7867_2.png) [@Kokora](https://discourse.julialang.org/u/Kokora)\
**Post date:** [December 9, 2022, 6:06pm UTC](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907/20 "2022-12-09T18:06:07Z")

</div>

@lrnv Thank you for your comment. I seem to still run into issues even when I use your snippet of code.  
I have another idea. What if I define the copula based on Y1 and a latent variable Y2st = logistic(mu2, 1). Also I supplement the likelihood with the Ycat (ordinal categorical data) using the OrderedLogistic likelihood. I tried to implement that below. But it does not seem to work either. I might be doing something terrible wrong here.

```julia

#-- Clayton Copula Model 180 Rotate
@model function Model3Clay180(Y1, Y2, X)

    n, K = size(X)
    ncateg = 4#maximum(Y2)
     
    # Priors
    cutpoints ~ Bijectors.ordered(filldist(Normal(0., 10.), ncateg - 1))
    θ ~ Normal(0.,1.)
    θ₁ ~ Normal(0.,1.)
    θ₂ ~ Normal(0.,1.)

    β1 ~ filldist(Normal(0., 1.), K)
    β2 ~ filldist(Normal(0., 1.), K)
     β ~ filldist(Normal(0., 1.), K)

    # Regression parameters
    μ1 = zeros(typeof(θ₁), n)
    μ2 = zeros(typeof(θ₂), n)
     ρ = zeros(typeof(θ), n)
     Y2st ~ filldist(Logistic(0, 1.), n)

    # Implement Shrinkage Priors
    for i in 1:n
        μ1[i] = θ₁ .+ sum(X[i,:].* β1)
        μ2[i] = θ₂ .+ sum(X[i,:].* β2)
        ρ[i] = exp(θ .+ sum(X[i,:].* β))
        Y2[i] ~ OrderedLogistic(Y2st[i]+μ2[i],cutpoints)

        #Now the copula is between Y1 (observed) and Y2st (latent)
        Turing.Turing.@addlogprob! loglikelihood(SklarDist(SurvivalCopula(ClaytonCopula(2,ρ[i]),(1,2)), (Normal(μ1[i], 1.), Logistic(μ2[i],1.))), [Y1[i],Y2st[i]+μ2[i]] )

    end
end
       

```

[Next page](https://discourse.julialang.org/t/issues-implementing-copula-regression-using-copula-and-turing/90907.md?page=2)