Turing: no method matching MvNormal(::Matrix{Float64}, ::Float64)

Specific Domains Statistics
1

Hello everybody,

as I want to switch from R to Julia for my Bayesian modelling, tonight I run my first models in Turing. However, I encountered some problems, principally when trying to model a random-intercept model.

The model is defined as follows:

@model my_regression(x, y, idx) = begin
    α ~ Normal(mean(y), 2.5 * std(y))
    β ~ Exponential(2)
    σ ~ Exponential(1 / std(y))

    n_gr = length(unique(idx))

    τ ~ truncated(Cauchy(0, 1), 0, Inf)     # group-level SDs intercepts
    αⱼ ~ filldist(Normal(0, τ), n_gr)           # group-level intercepts

    ŷ = α .+ αⱼ[idx] .+ x * β 
    y ~ MvNormal(ŷ, σ)
end;

However, when running the model I receive the error no method matching MvNormal(::Matrix{Float64}, ::Float64).

A similar model without the group-level intercepts runs without any problems and uses @. y ~ Normal(exp(α + β1 * x), σ) as a linear model. I suspect that a main reason for my problems could be due to a lack of understanding of the role of broadcasting, but I am really not sure.

Looking forward to your comments,

Tarotis

2

I could be wrong but i think MvNormal is expecting a vector and not a matrix as you appear to be sending it.

3

You are correct. ŷ should be a vector. I’m guessing x is a n by 1 array rather than a column vector.

4

What the others have pointed out is correct. There is a MatrixNormal, but that’s likely overkill here. This should be fine:

y ~ MvNormal(vec(ŷ), σ)

Edit: this assumes you also vectorize y.
Also, you’ll get better performance if you precompute constants like mean(y), std(y), and unique(idx) outside of the model. And this is faster:

ŷ = α .+ view(αⱼ, idx) .+ x .* β
2 Likes
5

Thank you all for the feedback and explicit explanations! It works great now. As the models will be quite complex in the end I fear that STAN/brms might be faster still, but I guess I will only find out once I have the corresponding models ready.