"sample from the prior" vs "sample from generative model'

[scheidan]

I'm trying to understand how to sample from the prior and how this contrasts to sampling from a "generative model'.

Do I interpret the example from the documentation correctly?

# define joint distribution p(s, m, x, y)
@model gdemo(x, y) = begin
    s ~ InverseGamma(2,3)
    m ~ Normal(0,sqrt(s))
    x ~ Normal(m, sqrt(s))
    y ~ Normal(m, sqrt(s))
    return x, y
end

# this samples from p(x,y)  
g_prior_sampler = gdemo()
g_prior_sampler()

# this samples p(x|y=-3.4)
c = sample(gdemo(missing, -3.4), NUTS(5000, 0.65))

Based on this I played around with the second example:

# define p(s, m, [x_1, ..., x_n])
@model gdemo2(x) = begin
    s ~ InverseGamma(2,3)
    m ~ Normal(0, sqrt(s))
    for i in eachindex(x)
        x[i] ~ Normal(m, sqrt(s))
    end
    return x
end

# sample from  p(s, m, x)
g_prior_sampler2 = gdemo2()
g_prior_sampler2()   # -> ERROR: MethodError: no methods matching keys(::Symbol)
# I guess the error is because the length of x in not defined.
# However passing the length as extra parameter does not resolve this.

# sample from p(s, m, [x_1, x_2, x_3])  (that's also sampling from the prior)
c = sample(gdemo2([missing, missing, missing]), NUTS(5000, 0.65))
# -> works :)

# sample from p(s, m, [x_1, x_2] | x_3=2.2)
c = sample(gdemo2([missing, missing, 2.2]), NUTS(5000, 0.65))
# -> ERROR: MethodError: no methods matching (::Colon)(::Int64, ::Base.OneTo{Int64})

The last case would be very useful in applications with missing data.

(probably related to #446.)

[cpfiffer]

You can fix the first case by passing in a default vector type for x, as below:

@model gdemo2(x=repeat(Union{Missing, Float64}[missing], 3)) = begin
    s ~ InverseGamma(2,3)
    m ~ Normal(0, sqrt(s))
    for i in eachindex(x)
        x[i] ~ Normal(m, sqrt(s))
    end
    return x
end

This allows Turing's compiler to infer both type and length for x.

As to your last example: I thought I remembered that we removed support (see #651) for vectors that are of type Union{Missing, Float64}, like what you have here: [missing, missing, 2.2]. There was a lot of conversation on this at #544, as well.

@mohamed82008 do you remember adding this back in? I can't seem to find it anywhere but I have a possibly false memory of you building this after #651.

[scheidan]

Thanks a lot! That works nicely.

More fundamentally, could someone confirm or correct if my probabilistic interpretations (see comments in the code example) are correct?

I'd be help to clarify this in the documentation.

[mohdibntarek]

No, the mixed missing case is not in there.

[cpfiffer]

More fundamentally, could someone confirm or correct if my probabilistic interpretations (see comments in the code example) are correct?

Yes, I believe that's correct -- but could someone else confirm that? @yebai @xukai92

I think it's worth adding those comments and clarifying this behavior to the documentation, I think they're a very helpful way to think about things. Would you mind if I put something similar to your comments in the documentation?

[scheidan]

I think that would be very helpful. For someone like me coming from statistics it is not exactly clear what sampling from the "generative model" means (also the ML community seems to use different
definitions: https://en.wikipedia.org/wiki/Generative_model). Also "sampling from the prior" could be interpreted more narrow than what the current implementation is capable of.

That's the summary of my understanding:

1. Sampling from a unconditional distribution. This could be the whole joint distribution p(a, b, c, d) or a marginal of it, e.g. p(a,b) = \int p(a, b, c, d) dcdb.

   This case is currently covered with

   prior_sampler = model()
   prior_sampler()

   if model defines a return value and your trick in case we have vectors.

   Alternatively, we can simply use sample by setting all arguments to missing.

2. Sampling from a conditional distribution. For example p(a, b | c=1.1, d=2.2), or if we assume d is a vector p(a, b, [d_2, d_4] | c=1.1, [d_1 = 3.1, d_3=3.3] ).

   For conditioning we need a sampler (HMC, ...), so it must be used via sample (with the current caveat that missings are not allowed in arrays).

Because both cases are covered by sample, it might be a good idea to remove the special "prior sampling" method.

Let me know if I can help you with the documentation.

[yebai]

More fundamentally, could someone confirm or correct if my probabilistic interpretations (see comments in the code example) are correct?

Yes, I believe that's correct -- but could someone else confirm that? @yebai @xukai92

@scheidan I can confirm the interpretation is correct. While "generative model" mostly refers to the "prior" in the ML community, this terminology a bit abused in my opinion. I think the unconditional / conditional distribution terminology is more self-explanatory and easier to work with.

Regarding "missing" data, we do want to support it in the future. It would be great to have some more concrete discussions on possible use cases first.

[xukai92]

I agree with Hong's comments.

Regarding

Regarding "missing" data, we do want to support it in the future. It would be great to have some more concrete discussions on possible use cases first.

For the example you provided

# sample from p(s, m, [x_1, x_2] | x_3=2.2)
c = sample(gdemo2([missing, missing, 2.2]), NUTS(5000, 0.65))
# -> ERROR: MethodError: no methods matching (::Colon)(::Int64, ::Base.OneTo{Int64})

Are you trying to impute the missing data via sampling from p([x_1, x_2] | s, m) during the sampling of p(s, m | x_3=2.2)?

[scheidan]

Are you trying to impute the missing data via sampling from p([x_1, x_2] | s, m) during the sampling of p(s, m | x_3=2.2)?

Yes exactly, that's just another way to look at it, because
p([x_1, x_2] | s, m) * p(s, m | x_3=2.2) = p(s, m, [x_1, x_2] | x_3=2.2).

As soon I have time I will provide some more explicit use-cases.

[xukai92]

Yes that's the same equation I was following. So the decision here is whether or not to provide a "default missing data imputation mechanism" as this? I guess the answer for now is no and we'd like to ask users to do this more explicitly after they obtain the chain.

[cpfiffer]

The guide is a hopefully a little clearer on this now. Does anyone have anything else they think I should put in there or anything that should be expanded on a bit more?

[scheidan]

Thanks, looks great!
In the section "Sampling from an Unconditional Distribution (The Prior)" it should be mentioned, that we must have a return statement of all variables we whish to sample.