# Issue with MvNormal and Turing
**URL:** https://discourse.julialang.org/t/issue-with-mvnormal-and-turing/54541
**Category:** Probabilistic Programming
**Created:** [February 3, 2021, 3:14pm UTC](https://discourse.julialang.org/t/issue-with-mvnormal-and-turing/54541 "2021-02-03T15:14:34Z")
**Posts on this page:** 9
**Page:** 1
### Author:  [@Carter](https://discourse.julialang.org/u/Carter)
#### Post date: [February 3, 2021, 3:14pm UTC](https://discourse.julialang.org/t/issue-with-mvnormal-and-turing/54541/1 "2021-02-03T15:14:34Z")
I’m trying to recreate the model from chapter 5.1 to estimate correlation from the book [Bayesian Cognitive Modeling: A Practical Course](https://cpb-us-e2.wpmucdn.com/faculty.sites.uci.edu/dist/0/180/files/2011/03/BB_Free.pdf) using Turing. However, I keep running into an issue with MvNormal where it does not like the covariance matrix because it says that it is not positive semi-definite and that Cholesky decomposition has failed or a LinearAlgebra.SingularException(2) error.
The code to reproduce these errors:
```julia
using Turing, Distributions, Random, Statistics
using Plots, StatsPlots, MCMCChains
dataset1 = [0.8 102; 1 98; 0.5 100; 0.9 105; 0.7 103; 0.4 110; 1.2 99; 1.4 87; 0.6 113; 1.1 89;
1.3 93]
dataset2 = [dataset1; dataset1]
println("Number of observations in dataset 1: ", size(dataset1, 1))
println("Number of observations in dataset 2: ", size(dataset2, 1))
#scatter(dataset1[:, 1], dataset1[:, 2], xlabel="Response Time (sec)", ylabel="IQ", title="Dataset 1", legend=false)
@model corr_model(x) = begin
#Prior
r ~ Uniform(-1, 1)
mu ~ MvNormal(zeros(2), 1/sqrt(0.001)) # performance tip from Turing tutorials
lambda = Vector{Real}(undef, 2)
lambda[1] ~ Gamma(0.001, 0.001)
lambda[2] ~ Gamma(0.001, 0.001)
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
model1 = corr_model(dataset1)
chain1 = sample(model1, NUTS(), 5000, progress=true)
density(chain1[:r], xlabel="Correlation", ylabel="Posterior Density", legend=false)
```
Any help would be appreciated.
---
### Author:  [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [February 3, 2021, 3:17pm UTC](https://discourse.julialang.org/t/issue-with-mvnormal-and-turing/54541/2 "2021-02-03T15:17:44Z")
Similar to [Bivariate normal in Turing - #5 by drbenvincent](https://discourse.julialang.org/t/bivariate-normal-in-turing/53446/5)
---
### Author:  [@drbenvincent](https://discourse.julialang.org/u/drbenvincent)
#### Post date: [February 3, 2021, 3:54pm UTC](https://discourse.julialang.org/t/issue-with-mvnormal-and-turing/54541/3 "2021-02-03T15:54:47Z")
I’m currently translating all the models from this book. Are you doing the same thing, or just focussing on that one example?
---
### Author:  [@Carter](https://discourse.julialang.org/u/Carter)
#### Post date: [February 3, 2021, 3:59pm UTC](https://discourse.julialang.org/t/issue-with-mvnormal-and-turing/54541/4 "2021-02-03T15:59:25Z")
I’m just focusing on that one example. Looking at your problem it seems that originally you were using std instead of your variances and when you fixed that it worked for you?
---
### Author:  [@drbenvincent](https://discourse.julialang.org/u/drbenvincent)
#### Post date: [February 3, 2021, 4:16pm UTC](https://discourse.julialang.org/t/issue-with-mvnormal-and-turing/54541/5 "2021-02-03T16:16:47Z")
Yep. The Normal distribution in JAGS is parameterised by mean and precision, but in most other approaches it’s by mean and std. So I’ve just been setting priors over std instead of exactly duplicating what’s in the text.
This is my full (and working) code. You can probably come up with more sensible priors than uniform over std but…
```julia
using Turing, StatsPlots
x = [0.8 102; 1 98; 0.5 100; 0.9 105; 0.7 103; 0.4 110; 1.2 99; 1.4 87; 0.6 113; 1.1 89; 1.3 93]
@model function pearson_correlation(x)
rows = size(x, 1)
# priors
μ₁ ~ Normal(0, 100)
μ₂ ~ Normal(0, 100)
σ₁ ~ Uniform(0, 100)
σ₂ ~ Uniform(0, 100)
r ~ Uniform(-0.99, 0.99)
Σ = [σ₁*σ₁ r*σ₁*σ₂;
r*σ₁*σ₂ σ₂*σ₂]
# likelihood
for i = 1:rows
x[i,:] ~ MvNormal([μ₁, μ₂], Σ)
end
end
chain = sample(pearson_correlation(x), NUTS(), 10000)
```
---
### Author:  [@Carter](https://discourse.julialang.org/u/Carter)
#### Post date: [February 3, 2021, 4:28pm UTC](https://discourse.julialang.org/t/issue-with-mvnormal-and-turing/54541/6 "2021-02-03T16:28:51Z")
Thanks for the link Mohamed. It’s the same problem set up in a different way. Unfortunately, I still can’t see why my model is causing the error in the covariance matrix.
---
### Author:  [@drbenvincent](https://discourse.julialang.org/u/drbenvincent)
#### Post date: [February 3, 2021, 4:50pm UTC](https://discourse.julialang.org/t/issue-with-mvnormal-and-turing/54541/7 "2021-02-03T16:50:41Z")
Perhaps stepwise convert the working code into what you have to see what bit causes the issue?
---
### Author:  [@jkbest2](https://discourse.julialang.org/u/jkbest2)
#### Post date: [February 3, 2021, 6:42pm UTC](https://discourse.julialang.org/t/issue-with-mvnormal-and-turing/54541/8 "2021-02-03T18:42:25Z")
I suspect the issue is the `Gamma` priors here. If you check `rand(Gamma(0.001, 0.001))`, you get a lot of very small or underflow-to-zero results. The `Distributions.jl` `Gamma` distribution is parameterized with shape and scale (i.e. the k and \theta parameterization listed on the [Wikipedia page](https://en.wikipedia.org/wiki/Gamma_distribution)). You probably want your `lambda` priors to be `Gamma(0.001, 1 / 0.001)` so that it has mean one.
---
### Author:  [@Carter](https://discourse.julialang.org/u/Carter)
#### Post date: [February 3, 2021, 7:29pm UTC](https://discourse.julialang.org/t/issue-with-mvnormal-and-turing/54541/9 "2021-02-03T19:29:07Z")
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.
```julia
@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
```