# Hyper Priors in Turing

**URL:** <https://discourse.julialang.org/t/hyper-priors-in-turing/20752>\
**Category:** Statistics\
**Tags:** package\
**Created:** [February 13, 2019, 4:00pm UTC](https://discourse.julialang.org/t/hyper-priors-in-turing/20752 "2019-02-13T16:00:53Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![krenova](https://avatars.discourse-cdn.com/v4/letter/k/53a042/32.png) [@krenova](https://discourse.julialang.org/u/krenova)\
**Post date:** [February 13, 2019, 4:00pm UTC](https://discourse.julialang.org/t/hyper-priors-in-turing/20752/1 "2019-02-13T16:00:54Z")

</div>

Hi,

I’m wondering if hyperpriors work on Turing. I’m playing around with Turing and trying to use it to estimate a Hidden Markov Model. In my HMM model, I’ve used a multivariate normal prior on parameters of my transition probability model with its mean and variance components respectively having a multivariate normal and inverse wishart prior (hyperpriors). With the hyper priors, I always get the error:  
`ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed.`

Its tough to explain the code and model and be concise so I’ve created a similar example on bayesian linear regression with similar prior and hyper priors to estimate the regression coefficients. This time however, I did not get the earlier error seen above but was instead returned a long error as follows:

`ERROR: MethodError: no method matching PDMats.PDMat{ForwardDiff.Dual{ForwardDiff.Tag{getfield(Turing.Core, Symbol("#f#24")){Turing.Core.VarReplay.VarInfo,Turing.Model{Tuple{:σ₂,:Sigma,:Mu,:beta_hat},Tuple{:y},getfield(Main, Symbol("###inner_function#31011#140")).......... and more...`

The code for the bayesian regression is as follows:

```julia
####################################################################
# Data Generation
####################################################################
df = DataFrame(
    x1 = rand(Normal(), 100), 
    x2 = rand(Normal(), 100)
);
beta = [5, 4, 5];
df[:y] = beta[1] .+ beta[2]*df[:x1] + beta[3]*df[:x2] + rand(Normal(), 100);

####################################################################
# Bayesian Inference with Turing
####################################################################
@model BayesReg(y, x1, x2) = begin
    
    n_obs = length(y)

    σ₂ ~ InverseGamma(0.001, 0.001)

    Sigma ~ InverseWishart(
        3 + 1, 
        Matrix{Float64}(I, 3, 3) 
    )
    Mu ~ MvNormal( 
        zeros(3),    
        Matrix{Float64}(I, 3, 3)   
    )

    beta_hat ~ MvNormal(Mu, Sigma)

    mu = beta_hat[1] .+ beta_hat[2]*x1 .+ beta_hat[3]*x2
    for i = 1:n_obs
        y[i] ~ Normal(mu[i], σ₂)
    end
end;

model = BayesReg(df[:y], df[:x1], df[:x2])
chain = sample(model, NUTS(1500, 200, 0.65));

```

Any ideas where I may have gone wrong or misunderstood something about the package…?

---

<div class="post-metadata">

**Author:** ![ElOceanografo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eloceanografo/32/624_2.png) [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)\
**Post date:** [February 13, 2019, 5:16pm UTC](https://discourse.julialang.org/t/hyper-priors-in-turing/20752/2 "2019-02-13T17:16:54Z")

</div>

I’m guessing the first error happens because the covariance matrix of your hyperprior is not quite symmetric due to floating point rounding error:

```julia
using Distributions, Random, LinearAlgebra
seed!(1)
A = randn(3,3);
A = A'A
B = randn(3,3)
B = B'B
C = A' * B * A
# Result looks symmetric...
# 3×3 Array{Float64,2}:
# 17.0807 -5.12167 -11.5349
# -5.12167 4.1515 5.32059
# -11.5349 5.32059 9.48419
# ...but isn't actually:
issymetric(C)
# false
MvNormal(zeros(3), C)
# ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed.

```

The solution is to wrap `C` in `Symmetric` to tell `MvNormal` it can ignore the floating-point errors and assume the matrix is in fact symmetrical:

```julia
MvNormal(zeros(3), Symmetric(C)) # works

```

This is a confusion that I’ve run into before, as have [other people](https://github.com/JuliaStats/Distributions.jl/issues/791). It would be great to get at least a more informative error message for it.

The second error is coming from the automatic differentiation done by the NUTS algorithm–if you try

```julia
chain = sample(model, MH(1500))

```

it runs. I can’t parse that error message either, but it looks like there’s some unhappy interaction going on between types from PDMats, ForwardDiff, and Turing.

---

<div class="post-metadata">

**Author:** ![krenova](https://avatars.discourse-cdn.com/v4/letter/k/53a042/32.png) [@krenova](https://discourse.julialang.org/u/krenova)\
**Post date:** [February 14, 2019, 1:09am UTC](https://discourse.julialang.org/t/hyper-priors-in-turing/20752/3 "2019-02-14T01:09:20Z")

</div>

Thanks @ElOceanografo, you are completely right about the issue with symmetry.

You are also completely right about the issue with the automatic differentiation (AD). I managed to resolve it by changing the AD to Turing.setadbackend(:reverse\_diff).

Can’t thank you enough, your answer was prompt and came with examples. Was about to give up on Turing!

---

<div class="post-metadata">

**Author:** ![ElOceanografo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/eloceanografo/32/624_2.png) [@ElOceanografo](https://discourse.julialang.org/u/ElOceanografo)\
**Post date:** [February 14, 2019, 1:26am UTC](https://discourse.julialang.org/t/hyper-priors-in-turing/20752/4 "2019-02-14T01:26:56Z")

</div>

No worries, glad you sorted it out!

Pinging @andreasnoack and @jrevels also…can you tell if the issue with ForwardDiff/PDMats here is something expected, or is it a bug?

---

<div class="post-metadata">

**Author:** ![mohamed82008](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mohamed82008/32/18171_2.png) [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)\
**Post date:** [February 14, 2019, 3:17am UTC](https://discourse.julialang.org/t/hyper-priors-in-turing/20752/5 "2019-02-14T03:17:59Z")

</div>

I found the problem, but while fixing it, I may have stumbled upon a bigger problem in Julia itself; it silently exits. I will open an issue in Julia.

Edit: [Julia 1.1 silently exits · Issue #31067 · JuliaLang/julia · GitHub](https://github.com/JuliaLang/julia/issues/31067)

---

<div class="post-metadata">

**Author:** ![Elmo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/elmo/32/17979_2.png) [@Elmo](https://discourse.julialang.org/u/Elmo)\
**Post date:** [February 17, 2021, 9:34pm UTC](https://discourse.julialang.org/t/hyper-priors-in-turing/20752/6 "2021-02-17T21:34:16Z")

</div>

Hi all,

Sorry to chime in here, but I am having the same issue. However, I’ve found that Turing only throws `ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed.` when the dimensionality of the covariance matrix increases above 3. For exampe, I tried @krenova’s MWE (3 dimensional Sigma) and everything works. However, when I slightly modify the code to include a 4th variable I get the PosDefException. Below is the modified code, including the `Symmetric(...)` fix suggested by @ElOceanografo.

```julia
df = DataFrame(x1 = rand(Normal(), 100), x2 = rand(Normal(), 100), x3 = rand(Normal(), 100))
beta = [5, 4, 5, 6]
df[!, :y] = beta[1] .+ beta[2]*df[!, :x1] + beta[3]*df[!, :x2] + beta[4]*df[!, :x3] + rand(Normal(), 100)

@model BayesReg(y, x1, x2, x3) = begin
    
    n_obs = length(y)

    σ₂ ~ InverseGamma(0.001, 0.001)

    Sigma ~ InverseWishart(
        4 + 1, 
        Matrix{Float64}(I, 4, 4) 
    )
    Mu ~ MvNormal( 
        zeros(4),    
        Matrix{Float64}(I, 4, 4)   
    )

    beta_hat ~ MvNormal(Mu, Symmetric(Sigma))

    mu = beta_hat[1] .+ beta_hat[2]*x1 .+ beta_hat[3]*x2 .+ beta_hat[4]*x3 
    for i = 1:n_obs
        y[i] ~ Normal(mu[i], σ₂)
    end
end

model = BayesReg(df[!, :y], df[!, :x1], df[!, :x2], df[!, :x3])
chain = sample(model, NUTS(0.65), 1000)

```

I’m using Turing v0.15.10. Does anybody know how to fix this? Any help would be greatly appreciated!

Ps.[The issue](https://github.com/JuliaStats/Distributions.jl/issues/791) mentioned by @ElOceanografo contains a suggestion that [GitHub - timholy/PositiveFactorizations.jl: Positive-definite "approximations" to matrices](https://github.com/timholy/PositiveFactorizations.jl) can be used to help with the apparent numerical issues, but I’m not sure how to actually do that within Turing…