Issue with MvNormal and Turing

I think I am one step closer to a solution and @jkbest2 has the right idea. For Gamma distributions, I believe Distributions.jl parameterizes the Gamma distribution with shape \alpha and scale \theta while WinBUGS uses shape \alpha and mean parameter \mu = \alpha\theta. So the Gamma(0.001, 0.001) in WinBUGS is Gamma(0.001, 1) or Gamma(0.001) in Julia.

But now I am often getting warnings for numerical errors when I run this model.

@model function corr_model(x, ::Type{T}=Float64) where T
    #Prior
    r ~ Uniform(-0.999, 0.999)
    # mu ~ MvNormal(zeros(2), 1/sqrt(0.001))
    mu = Vector{T}(undef, 2)
    mu[1] ~ Normal(0, 1/sqrt(0.001))
    mu[2] ~ Normal(0, 1/sqrt(0.001))
 
    lambda = Vector{T}(undef, 2)
    lambda[1] ~ Gamma(0.001, 1)
    lambda[2] ~ Gamma(0.001, 1)
    sigma = 1 ./ sqrt.(lambda)
    
    Σ = [1/lambda[1] r*sigma[1]*sigma[2]; r*sigma[1]*sigma[2] 1/lambda[2]]
    
    for i = 1:size(x, 1)
       x[i, :] ~ MvNormal(mu, Σ) 
    end
end