# ReactiveMP: How to run linear model with multiple predictors and an intercept

**URL:** https://discourse.julialang.org/t/reactivemp-how-to-run-linear-model-with-multiple-predictors-and-an-intercept/85290
**Category:** Probabilistic Programming
**Tags:** question, package, bayesian-inference, forneylab
**Created:** [August 4, 2022, 10:01am UTC](https://discourse.julialang.org/t/reactivemp-how-to-run-linear-model-with-multiple-predictors-and-an-intercept/85290 "2022-08-04T10:01:47Z")
**Posts on this page:** 1
**Showing post:** 6

<div class="post-metadata">

### Author: ![albertpod](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/albertpod/32/31336_2.png) [@albertpod](https://discourse.julialang.org/u/albertpod)
#### Post date: [August 5, 2022, 2:26pm UTC](https://discourse.julialang.org/t/reactivemp-how-to-run-linear-model-with-multiple-predictors-and-an-intercept/85290/6 "2022-08-05T14:26:54Z")

</div>

Hi @EvoArt!

Thanks for trying out `ReactiveMP.jl`.

You are trying to impose a prior on the diagonal elements of the covariance matrix of your likelihood function.  
We’ve just created a [PR](https://github.com/biaslab/ReactiveMP.jl/pull/184) that will support this setup (although you’d need to use precision parameter instead of variance and `MvNormalMeanPrecision` likelihood instead of `MvNormalMeanCovariance`).

To run your model with the current master branch of `ReactiveMP.jl` you can impose an `InverseWishart` distribution on the matrix of your likelihood function. I understand you don’t want to have additional correlations in your observations, so our solution is not 100% match your intentions:

```julia
using GraphPPL, Rocket, ReactiveMP, LinearAlgebra

n = 250
m = 100

@model [default_factorisation = MeanField()] function linear_regression(n,m)
    a ~ MvNormalMeanCovariance(zeros(m), diagm(ones(m)))
    b ~ NormalMeanVariance(0.0,1.0)
    W ~ InverseWishart(n+2, diageye(n))
    c ~ ones(n)*b
    x = datavar(Matrix{Float64})
    y = datavar(Vector{Float64})
    z ~ x*a+c
    y ~ MvNormalMeanCovariance(z, W)

end

results = inference(
    model = Model(linear_regression, n,m),
    data = (y = randn(n), x = randn(n,m)),
    initmarginals = (W = InverseWishart(n+2, diageye(n)), ),
    returnvars = (a = KeepLast(), b = KeepLast(), W = KeepLast()),
    free_energy = true,
    iterations = 10
);

using Plots
plot(results.free_energy) # check convergence

```

Note that you need to initialize marginals for the `W` (`initmarginals = (W = InverseWishart(n+2, diageye(n)), ),`) as some part of your computational graph will perform VMP. If you have further questions, we will be happy to help!

---

_[View the full topic](https://discourse.julialang.org/t/reactivemp-how-to-run-linear-model-with-multiple-predictors-and-an-intercept/85290)._