# Simple linear regression with ReactiveMP - RuleMethodError?

**URL:** https://discourse.julialang.org/t/simple-linear-regression-with-reactivemp-rulemethoderror/81526
**Category:** Probabilistic Programming
**Created:** [May 23, 2022, 6:29pm UTC](https://discourse.julialang.org/t/simple-linear-regression-with-reactivemp-rulemethoderror/81526 "2022-05-23T18:29:25Z")
**Posts on this page:** 1
**Showing post:** 2

<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: [June 8, 2022, 4:30pm UTC](https://discourse.julialang.org/t/simple-linear-regression-with-reactivemp-rulemethoderror/81526/2 "2022-06-08T16:30:24Z")

</div>

Hi @svilupp!

Thanks for trying out `ReactiveMP.jl`. We somehow missed your question on discourse; sorry for the long waiting reply.

Please note that assigning a Gamma prior to the precision parameter of Normal likelihood results in analytically intractable inference. In the context of message-passing, it means that you can’t execute belief propagation in your graph, which results in the first error: prefix `m_` stands for the _message_, while `q_` stands for _marginal_.

To circumvent this issue, we need to resort to, for example, variational inference or VMP in the context of factor graphs.

Now, there are a few ways to do that, but let me first point out that the only place where we need to factorize our model is just around `NormalMeanPrecision(a * time_index[i] + b, sigma)` factor. In particular, we need to use a mean-field factorization between the mean and precision, i.e.  
`a * time_index[i] + b` and `sigma`.

First and the easiest way to do that is to use the `default_factorisation` option when specifying the model:

```julia
@model [default_factorisation = MeanField()] function linreg(n)
    a ~ NormalMeanVariance(0.0, 10.0)
    b ~ NormalMeanVariance(0.0, 10.0)
    sigma ~ GammaShapeRate(1.0, 1.0)
    
    time_index = datavar(Float64, n)
    y = datavar(Float64, n)

    for i in 1:n
        y[i] ~ NormalMeanPrecision(a * time_index[i] + b, sigma)
    end
end

```

In this way, `ReactiveMP` will use mean-field wherever is possible, which in our case is just between `a * time_index[i] + b` and `sigma`.

Then the inference will follow smoothly:

```julia
results = inference(
    model = Model(linreg, length(time_index)),
    data = (y = y, time_index = time_index),
    initmessages = (b = vague(NormalMeanVariance),),
    initmarginals = (sigma = vague(GammaShapeRate),),
    returnvars = (a = KeepLast(), b = KeepLast(),sigma=KeepLast()),
    iterations = 20,
    showprogress = true,
)

```

About `initmarginals` and `initmessages`. Non-rigorously, you need to provide marginals when you resort to VMP. Likewise, you need to provide messages when the graph contains loops (technically there is a little more to that). Here, we need to initialize marginal for `sigma` and a message for either `a` or `b`.

Alternatively, you can create an auxiliary vector `aux` that will represent the means of your likelihood, i.e.:

```julia
@model function linreg(n)
    a ~ NormalMeanVariance(0.0, 10.0)
    b ~ NormalMeanVariance(0.0, 10.0)
    sigma ~ GammaShapeRate(1.0, 1.0)
    
    aux = randomvar(n)
    time_index = datavar(Float64, n)
    y = datavar(Float64, n)

    for i in 1:n
        aux[i] ~ a * time_index[i] + b
        y[i] ~ NormalMeanPrecision(aux[i], sigma)
    end
end

```

In this case, you would need to provide specific constraints on your posterior factorization:

```julia
constraints = @constraints begin 
    q(aux, sigma) = q(aux)q(sigma)
end

```

and feed them inside your inference function

```julia
results = inference(
    model = Model(linreg, length(time_index)),
    data = (y = y, time_index = time_index),
    constraints = constraints,
    initmessages = (b = vague(NormalMeanVariance),),
    initmarginals = (sigma = vague(GammaShapeRate),),
    returnvars = (a = KeepLast(), b = KeepLast(),sigma=KeepLast()),
    iterations = 20,
    showprogress = true,
)

```

This is a good extension of the Linear regression demo; please feel free to send PR with a new demo!

---

_[View the full topic](https://discourse.julialang.org/t/simple-linear-regression-with-reactivemp-rulemethoderror/81526)._