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.
3 Likes