Linear model with log-link on sigma (exp(sigma)) errors: are there plans to support it in the future?

[DominiqueMakowski]

I am trying to fit a simple linear model with some parameters expressed using a link-function (e.g., exp()) which is common to express priors in an unconstrained space.

using Random, RxInfer

y = randn(100)


RxInfer.@model function model_Gaussian(y)

    # Priors
    μ ~ RxInfer.NormalMeanVariance(0.3, 0.5)
    σ ~ RxInfer.NormalMeanVariance(log(0.2), 3)

    for i in eachindex(y)
        sigma = exp(σ)
        y[i] ~ RxInfer.NormalMeanVariance(μ, sigma)
    end
end

result = infer(
    model=model_Gaussian(),
    data=(y=y,),
)

Unfortunately, the above fails with:

ERROR: MethodError: no method matching exp(::GraphPPL.VariableRef{…})

In general, my question is about whether arbitrary link functions (such as that of StatsFuns) will be supported in the future?

Also, are there plans to make RxInfer work with a "standard" Distributions.Normal() rather than the bespoke NormalMeanVariance()?
Thanks for the clarifications!

[albertpod]

Hi @DominiqueMakowski. There's inherently no issue with RxInfer here, so I will transfer this issue into a discussion.
Given the non-conjugacy and nonlinearity of your model, you will likely have to use the CVI projections method. This is documented here.
Here's an example of how to solve this model with this method:

using RxInfer
using ExponentialFamilyProjection
using BayesBase
using Distributions

# Define the custom distribution first
struct TransformedNormalDistribution{H,T} <: ContinuousUnivariateDistribution
    h::H
    t::T
end

BayesBase.logpdf(dist::TransformedNormalDistribution, x) = logpdf(Normal(dist.h, exp(dist.t)), x)
BayesBase.insupport(dist::TransformedNormalDistribution, x) = true

@node TransformedNormalDistribution Stochastic [out, h, t]

# Define the model
@model function model_Gaussian(y)
    μ ~ Normal(mean=0.3, variance=0.5)
    σ ~ Normal(mean=log(0.2), variance=1.0)
    y .~ TransformedNormalDistribution(μ, σ)
end

# Define constraints
@constraints function non_conjugate_model_constraints()
    q(μ) :: ProjectedTo(NormalMeanVariance)
    q(σ) :: ProjectedTo(NormalMeanVariance)
    q(σ, μ) = q(σ)q(μ)
end

# Define initialization
@initialization function model_initialization()
    q(σ) = NormalMeanVariance(0.0, 1.0)
    # q(σ) = GammaShapeRate(1.0, 1.0) # also works
end

# Generate sample data
y = randn(100)

# Perform inference
result = infer(
    model = model_Gaussian(),
    data = (y = y,),
    constraints = non_conjugate_model_constraints(),
    initialization = model_initialization(),
    options = (rulefallback = NodeFunctionRuleFallback(),),
    showprogress = true,
    iterations = 10
)


println(result.posteriors[:μ][end])
println(result.posteriors[:σ][end])

I haven't checked the inference result, but this should give you a good starting point. There are other ways to run inference on your model, but I opted for the fastest one now.

As for support for Distributions.jl, I acknowledge it can sometimes be inconvenient. It's important to note that "primitive" nodes in GraphPPL are not the same as distributions from Distributions.jl. However, apart from model definition, RxInfer does support regular Distributions. We augment these with ExponentialFamily.jl, which provides some important features for RxInfer.jl's inference engine.

If you need any further clarification or have additional questions, feel free to ask.

[albertpod]

@Nimrais @bvdmitri Is it possible to avoid creating a fused node here and to combine the meta and constraints macros, so we don't change the appearance of the model?

[Nimrais]

@albertpod not now, as I am aware but we are working on it.

[bvdmitri]

@albertpod It is possible, but I needed to change the function, due to another issue with InverseGamma and variance.

@DominiqueMakowski here how you can do that:

transformation(σ) = inv(exp(σ)) # I use `inv` here because I mode, the precision instead, renders the inference more stable

RxInfer.@model function model_Gaussian(y)
    # Priors
    μ ~ RxInfer.NormalMeanVariance(0.3, 0.5)
    σ ~ RxInfer.NormalMeanVariance(log(0.2), 3)
    sigma := transformation(σ) # Use `:=`
    for i in eachindex(y)
        y[i] ~ RxInfer.NormalMeanPrecision(μ, sigma)
    end
end

Then we need to define approximation method for the transformation as well as the constraints and initial state to start variational iterations:

meta = @meta begin
    transformation(σ) -> CVIProjection()
end

constraints = @constraints begin 
    q(μ, σ, sigma) = MeanField()
    q(sigma) :: ProjectedTo(Gamma)
end

init = @initialization begin 
    q(μ) = NormalMeanVariance(0.3, 0.5)
    q(σ) = NormalMeanVariance(log(0.2), 3)
    q(sigma) = Gamma(1, 1)
end

all together:

result = infer(
    model=model_Gaussian(),
    data=(y=y,),
    constraints = constraints,
    meta = meta,
    initialization = init
)

runs without errors

Inference results:
  Posteriors       | available for (μ, σ, sigma)

and returns

result.posteriors[:σ] # NormalMeanVariance{Float64}(μ=0.4315748519638192, v=0.047652669642536616)

on my end.

[DominiqueMakowski]

Thanks for the thorough response. I feel like it's still a bit too complicated for me for now (this was mostly a simple toy example, and the solution contains quite a lot of extra steps that I would have no idea how to adapt to a more complex real-life scenario) - I need to learn more first before being able to fully utilize RxInfer ☺️

Maybe someday there will be a way to handle these steps somewhat automatically

It's important to note that "primitive" nodes in GraphPPL are not the same as distributions from Distributions.jl

That's a fair point, though I think many newcomer users might be confused seeing "Distributions.jl-like" syntax (as they might overlook the non-default named arguments) and expect things to just work, but I guess they'll have to read the docs :p

[albertpod]

@ismailsenoz seems like this is a simpler version of GCV node from here

[bvdmitri]

My question is whether arbitrary link functions (such as those in StatsFuns) will be supported in the future?

The original issue in your example arises because we use the := operator to deterministically map random variables, since we have limited control over the built-in Julia = operator. For example, you should write:

sigma := exp(...) # not sigma = exp(...) because the function involves random variables

I believe this is mentioned in the official GraphPPL documentation, but it is not covered in the RxInfer documentation. I’ll open a separate issue for that.

Can RxInfer work with a standard Distributions.Normal() rather than the custom NormalMeanVariance()?

No, not with a plain Normal. RxInfer supports multiple parameterizations of the Gaussian distribution, so you can use Normal(mean = ..., variance = ...), or for a shorter form, Normal(mean = ..., var = ...). The original Normal from Distributions actually has std as its second parameter and not variance, while MvNormal has covariance as its second parameters and not std. We believe its confusing so we opted to explicit parameters naming.

However, as @albertpod pointed out, you’ll still need to provide additional context for the inference engine. Even though the example might seem simple, the underlying inference is quite complex due to its non-conjugate structure. To handle this, you’ll need to write additional code. The best approach to solve the inference for this model is probably this example, which uses the exp function as desired (as Albert also already mentioned).

[wouterwln]

We have a note here about :=, but it might not be very clear. We might have to include a note in the documentation on a more prominent place that explains why = cannot be used to define factor nodes.