Error with MvNormal in Julia code for Bayesian GLM

[XinyueJia]

Hi everyone,

I'm still relatively new to Julia, Turing.jl, and Bayesian statistics, and I'm trying to implement a Bayesian GLM. However, when I try to run my code, I get an error message that says:

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

Here's the code that I'm running:

@model function bayesian_glm(y, G, P, L, A)
    # Priors for the coefficients
    β0 ~ Normal(0, 10)
    β1 ~ filldist(Normal(0, 1), 4)
    β2 ~ filldist(Normal(0, 1), 2, 2)
    β4 ~ Normal(0, 10)
    sigma ~ Exponential(1)

    μ = β0 .+ β1[G] .+ β2[L, P] .+ β4 .* A

    y ~ MvNormal(μ, sigma)
end
chain = sample(bayesian_glm(y, G, P, L, A), NUTS(), 1000, 4,  progress=false); 

The y variable is a continuous variable, G is a 4-factor categorical variable, and P, L, and A are 2-factor categorical variables. I'm trying to estimate the interactive effect between L and P.

I'm not sure what's causing this error, and I'm not familiar with the MvNormal function. Can anyone help me understand what's going wrong and how I can fix it?

Thanks in advance for any help you can provide!

[torfjelde]

MvNormal is provided by Distributions.jl and it was decided to disallow MvNormal(..., ::Real) because some users expected the 2nd argument (the sigma) to represent the standard deviation while others expected variance. Instead you can use sigma^2 * I where I, representing a uniform scale matrix, becomes available after doing using LinearAlgebra.

So, all in all

using LinearAlgebra

@model function bayesian_glm(y, G, P, L, A)
    # Priors for the coefficients
    β0 ~ Normal(0, 10)
    β1 ~ filldist(Normal(0, 1), 4)
    β2 ~ filldist(Normal(0, 1), 2, 2)
    β4 ~ Normal(0, 10)
    sigma ~ Exponential(1)

    μ = β0 .+ β1[G] .+ β2[L, P] .+ β4 .* A

    y ~ MvNormal(μ, sigma * I) # though if `sigma` is the variance, you should do `sigma^2 * I`
end
chain = sample(bayesian_glm(y, G, P, L, A), NUTS(), 1000, 4,  progress=false); 

[torfjelde]

Oooh my bad, I saw the scalar variance and immediately jumped on that. But then I'd recommend just doing

y ~ arraydist(Normal.(μ, sigma))

[XinyueJia]

Thank you! However, I got a weird error after running your codes: "DimensionMismatch: number of dimensions of x (1) must be greater than number of dimensions of d (2)"I have no idea what the x and d here refer to...

[torfjelde]

What shape is your y and A?

[devmotion]

It's for the same reason as before: y is a vector but you try to compare it with a distribution for matrices since μ is a matrix. I think the matrix-shape only arises from the interaction term β2[L, P]. Are you sure that it is correct and that you do not rather want getindex.((β2,), L, P)?

julia> A = rand(3, 2)
3×2 Matrix{Float64}:
 0.0694534  0.637922
 0.848685   0.205352
 0.800788   0.126694

julia> A[[1, 3], [2, 1]]
2×2 Matrix{Float64}:
 0.637922  0.0694534
 0.126694  0.800788

julia> getindex.((A,), [1, 3], [2, 1])
2-element Vector{Float64}:
 0.6379222569734143
 0.8007877016950632

[XinyueJia]

Oh cool! That makes a lot of sense! Thank you, David and Erlend, for your valuable help! (I'm new to Bayesian Stats and sorry for asking what may seem like a stupid question here... )

[XinyueJia]

Thanks, mister ChatGPT :D