Underestimated variance in the predictive message

[albertpod]

Following our discussion with @leostenzel and @bartvanerp, I've come across the following issue that occurs with nonmissing (observed) data points when the structured constraint is specified.

using RxInfer

data = (y = [1.0, 2.0, missing],) # (y = [missing, missing, missing],) works
@model function model_5(n)
    x = randomvar(n)
    y = datavar(Float64, n) where {allow_missing = true}

    z ~ NormalMeanPrecision(3, 100.0)

    γ ~ GammaShapeRate(1.0, 1.0)

    for i in 1:n
        x[i] ~ NormalMeanPrecision(z, 1.0)
        y[i] ~ NormalMeanPrecision(x[i], γ)
    end
end

@constraints function constraints_5()
    q(y, x, γ) = q(y, x)q(γ) 
end

result = inference(model = model_5(3), data = data, initmarginals = (γ=vague(GammaShapeRate), ), constraints = constraints_5(), iterations = 10,)

This results in

ERROR: MethodError: no method matching /(::PointMass{Float64}, ::Int64)

which makes sense given that we can't compute the mean_cov(q_out_μ).

To circumvent the error, we should use a different constraint:

@constraints function constraints_5()
    q(y, x, γ) = q(y)q(x)q(γ) 
end

which is a correct constraint but will result in underestimated variance in the prediction message (only γ will contribute).

EDIT: I wouldn't say it's necessarily a bug but it's something we need to think of.

[bvdmitri]

@albertpod This task has been added to the milestone for tracking and prioritization.

[albertpod]

Just dropping my thoughts on this issue.

When you specify {allow_missing = true}, RxInfer.jl starts treating the datavar as a random variable. Given the constraints, this effectively means that q_y_x contains the joint distribution between the point mass and a Gaussian. RxInfer.jl will attempt to trigger @rule NormalMeanPrecision(:τ, Marginalisation) (q_out_μ::Any,), which includes mean_cov(q_out_μ), and finally the function mean(f, itr) from Statistics.jl, where the division of this form occurs /(::PointMass{Float64}, ::Int64).

We can create a method for computing statistics of such joints. Alternatively, we can provide special treatment for datavars that contain both data and missing values and perform additional separation of the variables inside the infer function.

[bvdmitri]

I see. To properly fix this issue we probably should reimplement the way @rule macro works (which is by itself is a long standing issue). I can come up with a hot fix, but I don't think this would be ideal or work in all situations.