Dynamic models with index notation

[mhauru]

What is the statistical model that this should implement?

@model function f_indices()
    b ~ Bernoulli()
    x = Vector{Float64}(undef, b ? 1 : 2)
    if b
        x[1] ~ Normal(1.0)
    else
        x[1] ~ Normal(1.0)
        x[2] ~ Normal(2.0)
    end
end

Some options:

1. Nothing. That should be invalid.

2. A model for which the domain of x[2] is Union(Reals, {N/A}) and the joint is

p(b, x1, x2) = p(b) p(x1) p(x2 | b)

where

p(b) = Bernoulli
p(x1) = Normal(1.0)

and the distribution p(x2 | b) is such that its density function is

pdf(p(x2 | b == true))(y) = y == N/A ? 1.0 : 0.0
pdf(p(x2 | b == false))(y) = y == N/A ? 0.0 : pdf(Normal(2.0))(y)

3. A model like the above but where the domain of x[2] is Reals and

pdf(p(x2 | b))(y) = b ? 1.0 : pdf(Normal(2.0))(y)

This is what we currently do. I'm not convinced that this is a valid statistical model though.

Other similar models

@model function f_product()
    b ~ Bernoulli()
    if b
        dists = [Normal(1.0)]
    else
        dists = [Normal(1.0), Normal(2.0)]
    end
    x ~ product_distribution(dists)
end

This seems to currently implement option 2. from above.

@model function f_addinf()
    b ~ Bernoulli()
    x = Vector{Float64}(undef, b ? 1 : 2)
    if b
        x[1] ~ Normal(1.0)
    else
        if length(__varinfo__[@varname(x)]) == 1
            @addlogprob! -Inf
        end
        x[1] ~ Normal(1.0)
        x[2] ~ Normal(2.0)
    end
end

This, too, seems to implement something like option 2. from above, in a different manner. f_product errors if you evaluate it with a VarInfo where b is false and x only has one element, whereas f_addinf executes but rejects all such samples.

Sampling

Note that we don't know how to sample from f_product or from f_addinf, or more generally anything that implements option 2. from above. The only samplers that I know of that can handle both discrete and continuous parameters are MH and Gibbs. Gibbs won't ever work for this, because if x is in the acceptable range for some value of b then you can't get a valid sample from changing just b or changing just the range of x; You'd have to change both at once, and that is exactly what Gibbs doesn't do. MH could in principle work, but we would need to use a proposal distribution that allows going from x being length 1 to x being length 2. Such a jump would then have to be done at the same time as changing the value of b.

[penelopeysm]

The facts are that:

• the function returns a larger logp when b is true (x1 only) than when b is false (x1 and x2)
• MH acceptance probability is logp(x*)/logp(x)
• if x* contains b == true then it's more likely to be accepted
• thus the chain is biased towards b == true

There have been many arguments made about how this shouldn't be the case etc.

My point is that the way Turing is implemented now this is a logical and correct outcome. In short my view is that a Turing model is definedb by what happens when you run it, not an idealised mathematical equation that has a one-to-one correspondence with the model function. This is somewhat similar to operational and denotational semantics, so I'll use those terms in the following.

Even if one were to claim that Turing models do have denotational semantics (which is AFAIK not written / specified anywhere), it is clearly possible to construct Turing models that are meaningless:

@model function f()
    x ~ Normal()
    d = x < 0 ? Beta(2, 2) : Normal()
    x ~ d
end

This is mathematically bogus. To be fair, check_model will complain about this, but you can disable that and it will sample just fine. How do we attach meaning to the results of this? We can't understand it from a mathematical point of view (i.e. denotational), but we can understand it from the way it works (operational).

In the same way, the 'meaning' that I ascribe to the dynamic model is purely operational i.e. what happens when you pass it through a sampler), which is that it is biased towards fewer parameters.

At this point, I don't mind either of the following outcomes:

1. Status quo.
2. Declaring the dynamic model to be invalid (i.e. anything that does not have well-defined denotational semantics is out of scope).
3. Formalising a consistent denotational semantics of (dynamically-sized) Turing models and changing the implementation of MH to be in accord with that. This is where the word denotational breaks down a bit, I don't think it's necessary to formalise things in as much detail as the Wikipedia article suggests --- but one should have a logical, consistent, understanding of how Turing models should behave mathematically.

Note that describing a Turing model as a 'generative process' is also an operational definition. It's an alternative operational definition, which is implemented by calling rand(model). However, this is not consistent with what MH does and thus it's not applicable to sampling with MH. It's consistent with what Prior does.

On the other hand, claiming that MH and Prior should give the same result when there is no likelihood term is denotational in nature, because it expresses a property of a model regardless of what it looks like inside and without reference to how models are executed. Thus, if one were to assert this and change the behaviour of MH to work that way, then that would fall under (3).

And beyond that I have no real wish to argue for 1, 2, or 3.

[mhauru]

I would argue that we should have some denotational meaning attached to models because that is what makes PPLs powerful. A large part of the value of PPLs is in the abstraction between samplers and models. That you can define a model, a statistical model, and when you give it data to evaluate the likelihood with then you have implicitly also defined the posterior distribution of the model. And that any correct sampler for the posterior should yield values from that same distribution. If we lose this property then we are halfway back towards the situation where everyone who wants to do Bayesian inference using MCMC needs to understand how particular MCMC samplers work.

This abstraction will be imperfect. Maybe your Gibbs sampler can't be ergodic in your parameter space or something. But we should aim to keep the abstraction as strong as we can.

MH acceptance probability is logp(x*)/logp(x)

This is only true if the proposal distribution is symmetric. I'm a bit uncertain, but this might be relevant when thinking about proposal distributions between discrete sections of the domain of c2/x[2].

thus the chain is biased towards b == true

Doesn't this conclusion gloss over volumes of state space? Here's a model:

@model function f()
    x ~ Uniform()
    if x < 0.1
        @addlogprob!(0.0)
    else
        @addlogprob!(-4.4)
    end
end

The logpdf is higher in x < 0.1 but still the mean is about 0.1.