Standardising APIs between compiler and inference methods

[xukai92]

See https://github.com/TuringLang/Turing.jl/issues/634#issuecomment-471339521 for a summary.

[yebai]

@willtebbutt @xukai92
how about we always drop all observe statements in gdemo()(), ie, we always sample all variables from the prior, even for the following case raised by @willtebbutt (slightly modified)

@model demo(y) = begin
   5 ~ Normal(0,1)   # rv1
   y ~ Normal(0,1) . # rv2
end

gdemo()()  # sample both `rv1` and `rv2` from the prior, ingoring `rv1=5`
gdemo(2)() # sample from the posterior

This leads to consistent behaviour and should be intuitive too.

[xukai92]

I see. But what I mean by this issue is to actually have some concrete functions generated by evaluating gdemo() (I think @willtebbutt also proposed this some time ago).

Currently behaviour:

@model demo(y) = begin
   5 ~ Normal(0,1)   # rv1
   y ~ Normal(0,1) . # rv2
end
mf = gdemo()
mf() # => sampling from prior

But instead we could have

@model demo(y) = begin
   5 ~ Normal(0,1)   # rv1
   y ~ Normal(0,1) . # rv2
end
mf = gdemo(???) # this returns mf::UniqueRuntimeModelType
                # and create some function definitions:
logpdf(::UniqueRuntimeModelType, vi) = ... # (1) return logpdf by evaluating using vi
rand(::UniqueRuntimeModelType) = ... # (2) gdemo()() in your example
sample(::UniqueRuntimeModelType) = ... # (3) gdemo(2)() in your example

I guess it could also potentially give more space for performance because during sample we are explicitly using the specific version of logpdf and rand.

[willtebbutt]

@yebai your proposal has the issue that the behaviour of the following programme becomes ambiguous:

@model demo() = begin
    5 ~ Normal(0, 1)
end
demo() # what does this do?

This makes me nervous. On the one hand, this programme is weird, on the other it's a valid Turing programme, so we should ensure that the behaviour is well-defined. My proposal below will mostly resolve this.

@xukai92 I completely agree regarding the format of logpdf. Regarding rand and sample: certainly we could define one to sample from the prior and the other the posterior, but this wouldn't be particularly intuitive (there's nothing in the names to suggest the proposed behaviour). Instead, I propose the following:

1. Don't define sample, just use rand and logpdf (or maybe logp, since logpdf isn't defined for all models supported by Turing?).
2. The current proposal regarding gdemo(missing) to indicate that the prior should be used for y seems sensible to me.
3. We implement gdemo(Prior()) to indicate that all random variables should be sampled from the prior. In this case, both rv1 and rv2 would be sample from the prior.
4. rand(mf) produces a VarInfo. rand(mf, N) produces a Vector{<:VarInfo} of length N.
5. logpdf(mf, vi) works as proposed by @xukai92.

[mohamed82008]

One the one hand, this programme is weird, on the other it's a valid Turing programme, so we should ensure that the behaviour is well-defined.

We can give an error if the model has no variables.

UniqueRuntimeModelType

This can be just const PriorModel{pvars, F} = Model{pvars, Tuple{}, F, NamedTuple{(), Tuple{}}} which means no data, all parameters. For this case, no special sampling alg is required. More generally however, we can make use of the existing Sampler type, which we can add model (and vi) to. Then we can simply do: rand(Sampler(model, alg)) or for models without parameters, it can be rand(Sampler(model)).

Btw, after the new refactor, if you call gdemo(), it will give you a PriorModel{pvars, F}, with pvars as Tuple{:y}.

[willtebbutt]

We can give an error if the model has no variables.

Agreed, but this programme definitely has a random variable. It's anonymous, but it's definitely there, and if we define sample-from-prior functionality, then presumably something should happen.

Btw, after the new refactor, if you call gdemo(), it will give you a PriorModel{pvars, F}, with pvars as Tuple{:y}.

This is probably fine for now, but I'm not sure it's appropriate long-term for the same reasons as above.

This can be just const PriorModel{pvars, F} = Model{pvars, Tuple{}, F, NamedTuple{(), Tuple{}}} which means no data, all parameters. For this case, no special sampling alg is required. More generally however, we can make use of the existing Sampler type, which we can add model (and vi) to. Then we can simply do: rand(Sampler(model, alg)) or for models without parameters, it can be rand(Sampler(model)).

I like this. We would have to think carefully about how the rand(Sampler(model, alg)) would work for MCMC though. We should also probably implement a few convenience functions as well i.e. rand(model) = rand(Sampler(model)) and similar.

[mohamed82008]

It's anonymous, but it's definitely there, and if we define sample-from-prior functionality, then presumably something should happen.

Is there a practical application of sampling this anonymous random variable, or is it just for theoretical appeal? Note that implementation-wise, there is nothing random about a model with no explicit random variables; every time you run the model, you will get exactly the same thing, hence my difficulty digesting this anonymous random variable concept. I can see how 5 plays the role of data. So perhaps more generally this falls under likelihood models, where the only sensible thing to do with it, that I can think of, is to calculate the likelihood.

We should also probably implement a few convenience functions as well i.e. rand(model) = rand(Sampler(model)) and similar.

Sounds good.

[mohamed82008]

Note that currently, we just ignore any numbers or arrays on the LHS of a ~. If you want, we can assign an anonymous dvar to it instead and assign the value on the LHS of the ~ to this anonymous dvar. This would probably be more theoretically accurate. And it allows us to view the model with 5 ~ Normal(0, 1) only as a LikelihoodModel with 1 dvar.

Edit: by ignore I mean that we don't treat it as an explicit variable, but it still affects the logp calculation.

[willtebbutt]

Is there a practical application of sampling this anonymous random variable, or is it just for theoretical appeal?

I'm just looking for consistency. I don't have a particular application in mind :)

Note that implementation-wise, there is nothing random about a model with no explicit random variables; every time you run the model, you will get exactly the same thing.

Not if you run the model in sample-from-prior mode.

If you want, we can assign an anonymous dvar to it instead and assign the value on the LHS of the ~ to this anonymous dvar. This would probably be more theoretically accurate. And it allows us to view the model with 5 ~ Normal(0, 1) only as a LikelihoodModel with 1 dvar.

This general kind of thing would be a good change. If the model is constructed in the regular way, it should be a dvar as you suggest, but if the model is run in sample-from-prior mode, then it should presumably magically become a dvar.

[mohamed82008]

Not if you run the model in sample-from-prior mode.

@model demo() = begin
   5 ~ Normal(0,1)
end
model = demo()
logp = model(Turing.VarInfo(), Turing.SampleFromPrior())
all(isequal(logp), (model(Turing.VarInfo(), Turing.SampleFromPrior()) for i in 1:100))

vi is empty in every case, with a logp only.

[mohamed82008]

then it should presumably magically become a dvar.

This should be totally possible when constructing the Sampler.

[willtebbutt]

Sorry, I should have been clearer. I was referring to what I think should be going on rather than what actually happens at the minute. You're totally correct regarding the current behaviour. It's my feeling that if we implement a sampler-from-prior mode, then we should take these anonymous random variables more seriously and treat them more like named random-variables.

[mohamed82008]

then we should take these anonymous random variables more seriously and treat them more like named random-variables.

I don't foresee a problem handling these. I can try a proof-of-concept after the VarInfo refactor. Perhaps @trappmartin can also help with this after he reads the new compiler code.

[trappmartin]

Thanks for cc‘ing me. I’m pretty much in favour of the vi = rand(Sampler(model)) and the logpdf(model, vi) concept. Though, I believe it would be good if we keep the sample function. It’s semantically more clear to me what this does and that I can expect it to generate samples using a MCMC sampler. I would not expect it to draw from the prior as rand.

However, I’m not sure how inference with VI would fit into those syntaxes to stay consistent without having too many entry points for the user.

[willtebbutt]

However, I’m not sure how inference with VI would fit into those syntaxes to stay consistent without having too many entry points for the user.

My thinking here is that VI works fine with the rand(sampler(model, alg), N) design, where the alg parameter would be a representation of the approximate posterior distribution. Semantically this reads as: draw N samples from the posterior distribution over the random variables in model (each of which should be a VarInfo?) using alg. This works for both VI and MCMC, so provides a single point of entry for both kinds of samplers.

Though, I believe it would be good if we keep the sample function. It’s semantically more clear to me what this does and that I can expect it to generate samples using a MCMC sampler.

Again, my feeling here is that we wouldn't need the sample function, as the creation of the sampler, e.g. rand(sampler(model, HMC(...)), N), makes it quite clear that we're drawing N samples from model using the specified HMC algorithm. Although I could live with the sample function being an alias for rand, I suspect that rand is more consistent with existing stuff.

I’m pretty much in favour of the vi = rand(Sampler(model)) and the logpdf(model, vi) concept.

As am I :)

A minor stylistic thing, @yebai what would your thoughts be on making the "correct" thing to do be sampler(model, alg) rather than Sampler(model, alg)? (At least to me) the latter says make me an object of type Sampler whereas the former says make me a thing that can generate samples from model using alg, which may or may not be an object of type Sampler. The distinction is important because the Sampler option commits us to always producing an object of type Sampler if we want to produce samples, whereas the former allows us to spit out whatever object we like provided that it can be used to generate samples via rand.

[xukai92]

Reading through your explanation of sample v.s. Sampler, I feel sample is better, although in most of the cases (at this moment), sample will return a Sampler to us. Also I think using a function rather than a type conversion syntax is a more standard Julia style.