Hi there,

I’ve seen the recent development for the `LKJCholesky`

implementations within `Distributions.jl`

and am attempting to put it to use within `Turing.jl`

. I understand that `Turing.jl`

doesn’t yet support *Cholesky Variate* distributions. However, @sethaxen gave a few tips within the `Turing`

Slack channel on how to use it in `Turing.jl`

anyway:

Okay I underestimated how long a proper write-up will take, so here’s the short version. It uses TransformVariables instead of Bijectors, since that simplifies the code and lets us use only API functions.Instead of doing

`R ~ LKJ(n, eta)`

import TransformVariables and do

`trans = CorrCholeskyFactor(n) R_tilde ~ filldist(Flat(), dimension(trans)) R_U, logdetJ = transform_and_logjac(trans, R_tilde) Turing.@addlogprob! logpdf(LKJCholesky(n, eta), Cholesky(R_U)) + logdetJ`

When using the correlation matrix, instead of doing

`Σ = Diagonal(σ) * R * Diagonal(σ) y ~ MvNormal(μ, Σ)`

import PDMats and do

`Σ_L = LowerTriangular(collect(σ .* R_U')) # collect is necessary for ReverseDiff for some reason Σ = PDMat(Cholesky(Σ_L)) y ~ MvNormal(μ, Σ)`

Given this info I put together an example modeled after this blog post on Using a LKJ prior in Stan. It has simulated data but I am seeing `SingularException`

errors with this example. It doesn’t always error and sometimes samples really well.

Simulated data:

```
using Distributions, Turing, LinearAlgebra, Random
Random.seed!(121)
# generate data
sigma = [1,2,3]
Omega = [1 0.3 0.2;
0.3 1 0.1;
0.2 0.1 1]
Sigma = Diagonal(sigma) * Omega * Diagonal(sigma)
N = 100
J = 3
y = rand(MvNormal(Sigma), N)'
```

The model in `Turing.jl`

:

```
using TransformVariables, PDMats
# model
@model function correlation_cholesky(J, N, y)
sigma ~ filldist(truncated(Cauchy(0., 5.); lower = 0), J)
# LKJ work around
n, eta = J, 1.0
trans = CorrCholeskyFactor(n)
R_tilde ~ filldist(Flat(), dimension(trans))
R_U, logdetJ = transform_and_logjac(trans, R_tilde)
F = Cholesky(R_U)
Turing.@addlogprob! logpdf(LKJCholesky(n, eta), F) + logdetJ
Σ_L = LowerTriangular(collect((sigma .+ 1e-6) .* R_U'))
Sigma = PDMat(Cholesky(Σ_L))
for i in 1:N
y[i,:] ~ MvNormal(Sigma)
end
return (; Omega = Matrix(F))
end
```

Then sample the model:

```
mod = correlation_cholesky(J, N, y)
chain_chol = sample(mod, NUTS(0.65), 2000)
```

It’s here that the error commonly occurs. Any info would be appreciated.