Calling models on non-arguments

[cscherrer]

Say we have a model

m = @model μ begin
    x ~ Normal(0,1)
    y ~ Normal(x,1)
    z ~ Normal(y,1)
end

Currently we can do

julia> rand(m(μ=0))
(μ = 0, x = 0.7665341287079652, y = 2.641313887511112, z = 1.4301554808799748)

This works as expected, but a user might be surprised by

julia> rand(m(μ=0, y=10))
(μ = 0, x = 0.47837512610146454, y = 0.667663356964404, z = 1.388557652149238)

This happens because of the current rules for scoping. y is not an argument, so the model considers the first y to be outside the scope of the model.

QUESTION: Is this the semantics we want?

---

Note: You can still easily pass y into the model, like this:

julia> Do(m, μ=0,y=10) 
Joint Distribution
    Bound arguments: [μ, y]
    Variables: [z, x]

@model (μ, y) begin
        z ~ Normal(y, 1)
        x ~ Normal(μ, 1)
    end


julia> Do(m, μ=0,y=10) |> rand
(μ = 0, y = 10, z = 11.053045674424801, x = -1.6135437791883571)

So the concern isn't that this is hard to do (no pun intended). The point is more to check that the we have consistent reasonable semantics.

[millerjoey]

m = @model μ begin
  x ~ Normal(0,1)
  y ~ Normal(x,1)
  z ~ Normal(y,1)
end

Should this be x ~ Normal(μ,1)?

---

I think the most consistent thing to do for now would be to error and suggest a transformation that moves parameters to the args like predictive or Do.

The other option that's a little more afield is to do the transformations on-demand. In these cases though the user might need to distinguish between Doing and conditioning. Should MCMC be done too to sample the conditionals? Symbolic math if possible? Or just sample the aboves from predictive?

If that were possible, I might also expect that e.g. logpdf(m, (μ=1, x=0)) would be allowed and transform to something equivalent to logpdf(prior(m, :x), (μ=1, x=0)).

Should this all be done under the hood? It would be pretty interesting, but it's hard to know when to stop :P

[cscherrer]

Should this be x ~ Normal(μ,1)?

Good catch, yes that's what I meant.

The other option that's a little more afield is to do the transformations on-demand. In these cases though the user might need to distinguish between Doing and conditioning. Should MCMC be done too to sample the conditionals? Symbolic math if possible? Or just sample the aboves from predictive?

There's a core idea I think we need to preserve here: a Soss model is a joint distribution over its parameters, conditional on its arguments. This lets us keep model transformations and forward sampling very fast and simple. When we need to sample from some parameter conditional on observing something downstream, there are lots of options, and things get much more complicated.

So the guideline (at least for now) is that the core of Soss should be able to do things that are very mechanical, and where "the way to do it" is obvious and can be fast. Things with more choices (like inference) need to be outside the core.

Your example logpdf(m, (μ=1, x=0)) is really interesting here, because we can get from μ to x without any intermediate stochastic computations. I mean, there's no need to marginalize over anything. I'll need to think some more about the best way to do this. It's not hard, just a matter of the right design choice.

[millerjoey]

Yeah, I like the KIS approach in general. We can go a little further with "obvious" extensions if they're unambiguous and consistent.

Regarding calling models on non-arguments, here's another example to note, suggesting the idea of an argument might not be central:

julia> m = @model begin
       x ~ Normal(μ)
       y ~ Normal(x)
       end;
julia> rand(m(μ=4)) # this works as if μ is an arg
(x = 5.500516432453989, y = 4.208745378540177)

julia> m.args # but...
0-element Array{Symbol,1}

Likewise if m(x=0) would mean Do(x=0) (which is a nice idea IMO)also suggests the concept of a Model's args needs to be thought over.

[cscherrer]

Just to note, the current behavior has caused some confusion, e.g. @baggepinnen 's comment here.

[cscherrer]

Oh, and also...

julia> using Soss

julia> μ = 10.0
10.0

julia> m = @model begin
           x ~ Normal(μ)
           y ~ Normal(x)
       end;

julia> rand(m())
(x = 10.151203386020851, y = 10.506361236271212)

I think this is the right behavior, since otherwise things like Normal would either have to be passed as arguments or considered "special".