Simple model with a joint factor

[keorn]

I would like to specify a simple model, where I use a multivariate normal for one of the factors. The structure is:

Normal <-> MvNormal <-> Normal

However none of the below specifications succeed:

@model function model(a)
    a ~ Normal(mean = 1., var = 5.)
    b ~ Normal(mean = 2., var = 5.)
    
    [a, b] ~ MvNormal(mean = zeros(2), covariance = diageye(2))
end

MethodError: no method matching generate_get_or_create(::Expr, ::Nothing, ::Expr)

@model function model(ab)
    ab ~ MvNormal(mean = zeros(2), covariance = diageye(2))
    
    ab[1] ~ Normal(mean = 1., var = 5.)
    ab[2] ~ Normal(mean = 2., var = 5.)
end

Variable ab is already an individual variable in the model

@model function model(a)
    ab ~ MvNormal(mean = zeros(2), covariance = diageye(2))
    
    a ~ Normal(mean = 1., var = 5.)
    b ~ Normal(mean = 2., var = 5.)
    
    a ~ ab[1]
    b ~ ab[2]
end

Half-edge has been found: ab_1. To terminate half-edges 'Uninformative' node can be used.

Do you have any pointers to similar examples or general tips about the approach?

[albertpod]

Hi @keorn, could you please write the probabilistic generative model? I fail to see what are your observations here.

You could write something like this, but not sure that's what you mean.

using RxInfer

@model function ivy_model(y)
    a ~ Normal(mean = 1., var = 5.)
    b ~ Normal(mean = 2., var = 5.)
    
    c1 = [1, 0]
    c2 = [0, 1]
    y ~ MvNormal(mean = c1 * a + c2 * b, covariance = diageye(2))
end

result = infer(model=ivy_model(), data=(y=[10, 200],), free_energy=true)

[albertpod]

I think it's possible. Tagging folks for this @bvdmitri, @wouterwln, @bartvanerp

[bvdmitri]

@keorn thanks for trying out RxInfer! so of my thoughts, sorry for the rumble

I just specify the prior for now, no observations yet.

RxInfer requires at least some observations — it won’t work if none are provided.

ab ~ MvNormal(mean = zeros(2), covariance = diageye(2))

ab[1] ~ Normal(mean = 1., var = 5.)
ab[2] ~ Normal(mean = 2., var = 5.)

This is not possible in the current version of RxInfer. We’ve had quite a bit of debate about it, but there’s no clear consensus on how this should work in principle for arbitrary distributions (let alone how it could be implemented). For Gaussians, it might be possible to hack something together, but Albert’s proposal above is equivalent and should already work in the current version.

The key issue is that

ab ~ some_arbitrary_distribution(some_parameter = 1.0)

creates a single variable node ab in the graph. It is not a vector or a collection — it’s just one variable. Indexing like ab[1] only makes sense for collections of variables, not for a single variable node. It differs from other PPLs where ~ means "sample a value from RHS distribution", where you obtain an actual sample of size 2 which you can of course index as ab[1]. In RxInfer ~ really means "create a connection in the factor graph between LHS and RHS".

Since ab can be connected to arbitrary factor nodes that send arbitrary messages, there’s no way (at model construction time) to know its dimensionality, or even whether it’s univariate, multivariate, or something else. That information only becomes available at run time, but the factor graph must already be fully constructed at compile time.

I think the best shot is the example provided from Albert (works only for Gaussians)

the @model syntax supports only directed model specification

Not quite — our models are explicitly non-directed factor graphs. We don’t even store direction information. Or did you mean something else here?

[keorn]

RxInfer requires at least some observations — it won’t work if none are provided.

That does not seem to be the case, for an even simpler model Normal <-> Normal I am able to find the prior:

@model function model(a)
    a ~ Normal(mean = 1., var = 5.)
    a ~ Normal(mean = 2., var = 5.)
end
result = infer(
    model = model(),
	predictvars = (a = KeepLast(),)
)
mean(result.predictions[:a])

gives me 1.5 as expected.

What I meant with directionality is - it is not easy to specify the factors as joint distributions over the connected variables, one has to pick all variables except one and then define a distribution on that remaining variable.

Thanks for the clarification, I understand the challenges, will try things out using the methodology suggested by @albertpod .

[bvdmitri]

That does not seem to be the case, for an even simpler model Normal <-> Normal I am able to find the prior:

Ah, well, here you technically have an observation a (an argument to the model). Since you used it in the predictvars it will automatically pass missing as an observation to a. You can confirm that by doing

@model function model(a)
    a ~ Normal(mean = 1., var = 5.)
    a ~ Normal(mean = 2., var = 5.)
end
result = infer(
    model = model(), 
    data = (a = missing, )
)
mean(result.predictions[:a])

which should yield an identical result.

I'm not sure if we document this exact behavior well enough, but its definitely by design. If this is not present in the documentation we should definitely add it.

[keorn]

Right, so missing observation data is the same as no observation, makes sense.