is it or should it be possible to pass a distribution as data argument?

[skoghoern]

my attempt to use a distribution as data sample failed, when trying to build a hierarchical NormalMixture model.
I understand that the data variable should be a delta node. but currently it doesnt allow to pass a distribution. is this a misunderstanding on my side or something to implement?

@model function data_as_dist_example(x)
    z ~ NormalMeanVariance(0,1)
    x ~ NormalMeanVariance(z, 1)
end

result = infer(
    model = data_as_dist_example(),
    data = (x = NormalMeanVariance(1,1),),
    predictvars = (x = KeepEach(),)
)

#### error message ####
error: Stacktrace:       
[1] mean(fn:typeof(identity), distribution: :PointMass{MvNormalMeanPrecision{Float64, Vector{Float64}, Matrix{Float64}}})         @ BayesBase ~\.julia\packages\BayesBase\THarC\src\densities\pointmass.jl:29       
[2] mean(itr:PointMass{MvNormalMeanPrecision{Float64, Vector{Float64}, Matrix{Float64}}}) @ Statistics ~\.julia\packages\Statistics\gbcbG\src\Statistics.jl:44

#### debug info ####
@rule NormalMeanVariance (:μ, Marginalisation) (("marginal_out", "marginal_v")) - ReactiveMP.jl/src/message.jl
q_out: PointMass{NormalMeanVariance{Float64}}(NormalMeanVariance{Float64}(μ=1.0, v=1.0))
q_v: PointMass{Int64}(1)

Note: I added predictvars = (x = KeepEach(),) since otherwise we run into an error during predictvars and data setup in batch.jl (inference_check_dataismissing(d) = (ismissing(d) || any(ismissing, d)) ⇒ any(ismissing, NormalMeanVar) is not defined, since it needs an implementation of Base.iterate())

[meditans]

If I understand correctly, you are trying to pass an observation for the random variable x, that you want to be sampled from the NormalMeanVariance(1,1) distribution. In that case, I think you should do:

using RxInfer

@model function data_as_dist_example(x)
    z ~ NormalMeanVariance(0, 1)
    x ~ NormalMeanVariance(z, 1)
end

observation = rand(NormalMeanVariance(1, 1))

result = infer(
    model = data_as_dist_example(),
    data = (x = observation,),
)

The data you have to pass is just a datapoint. You're right in saying that in the inference process that datapoint will be converted in the factor $\delta(x - \text{observation})$, but data just takes the pure observation. You also don't need the predictvars option.

[skoghoern]

Hi Carlo,
thanks a lot for your answer. you made me rethink and reflect my question. i guess my initial question was actually quite confusing and not really making sense since i mixed up a message of a variable vs. a variable node. the solution you proposed is pretty much sampling from the distribution and using that as input data - thats a very neat idea, but not directly what i was looking for. my goal was not to use a pointmass but rather the combination of mean with variance to cluster several distributions.
i now understood i could simply pass it as a message by introducing another factor node, where we define the parameters with the "data":

@model function data_as_dist_example(m,v)
    z ~ NormalMeanVariance(0,1)
    x ~ NormalMeanVariance(z, 1)
    x ~ NormalMeanVariance(m,v)
end

result = infer(
    model = data_as_dist_example(),
    data = (m = 2,v=2) )

i actually found the solution by sketching out a factor graph, since in code it seems quite confusing to me. and i think my main help was the description from dimitri in #515 explaining that a the variables "defined" with the ~ arent really samples in rxinfer but rather "creating a connection between variable node and factor node".
quoted:

[x~dist(param)] creates a single variable node x in the graph. It is not a vector or a collection — it’s just one variable. Indexing like x[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".

[skoghoern]

and in case you or others are interested, thats a working version of the "hierachical normal mixture model" i was actually trying to implement:

#### working
@model function hierarchical_normal_mixture_two_levels(x1)
    local m1
    local w1
    local m2
    local w2

    p1 ~ Beta(1.0, 1.0)
    p2 ~ Beta(1, 10)
    for j in 1:2 # nr_components_level2
        m2[j] ~ NormalMeanVariance(0,1)
        w2[j] ~ Gamma(shape = 1, scale = 1)
    end
    for i in 1:2 # nr_components_level1
        s2[i] ~ Bernoulli(p2)
        x2[i] ~ NormalMixture(switch=s2, m=m2, p=w2)

        m1[i] ~ NormalMeanVariance(0,1)
        w1[i] ~ Gamma(shape = 1, scale = 1)
        x2[i] ~ NormalMeanPrecision(m1[i], w1[i])
    end
    s1 ~ Bernoulli(p1)
    x1 ~ NormalMixture(switch=s1, m=m1, p=w1)
end

result = infer(
    model = hierarchical_normal_mixture_two_levels(),
    data = (x1 = 1,),
    constraints = MeanField(),
    initialization = @initialization(begin  
    q(p2) = Beta(1,1)
    q(m2) = NormalMeanVariance(0,1)
    q(w2) = GammaShapeScale(1, 1)

    q(p1) = Beta(1,1)
    q(m1) = NormalMeanVariance(0,1)
    q(w1) = GammaShapeScale(1, 1)
    q(x2) = NormalMeanPrecision(0,1)
                     end),
)