Add Slab-and-Spike prior

[trappmartin]

It might be a good idea to add a slab-and-spike prior into Turing.

Currently, user can only do the following:

s ~ Bernoulli(p)
z ~ Normal(0, v)
b ~ s*z

In this case the user would have to use SMC, PG or a combination of PG and HMC/NUTS. A custom distribution would allow users to use HMC only, which would be a good feature. Something like:

b ~ SlabAndSpike(p, v)

I should have a look how this is handled in Stan link.

[richardreeve]

I could easily be wrong here as I've never used Stan, but as I understand it, they don't allow discrete distributions, so they can't have Dirac deltas making distributions discontinuous. As a result, they have a hack to avoid having a true spike in their slab-and-spike - I think they use a product of two or more normals(?) - product_normal() - that approximates a spike as the scale parameters tend to infinity... so it's not a real slab-and-spike either.

Although it would obviously be great to have a true slab-and-spike implemented directly in Turing.jl, my plan is to work on expanding the functionality of Distributions.jl to allow discontinuous distributions in general, and to help with the work to allow discrete distributions have support over any countable set rather than just integers - see JuliaStats/Distributions.jl#925 for both thoughts at the moment. Then the broader question would be how to allow these to be used in Turing, and in particular from the Distributions end, what information would we have to provide you and how for you to be able to automatically handle such distributions?

[trappmartin]

This is correct. I'll have to check what the product normal actually does.

A real slab-and-spike prior will anyhow only be usable if the user is using any particle-based algorithm. Other methods involve gradients, while SMC and co. only need to be able to sample from the prior using rand and to evaluate a samples likelihood using logpdf. There are no additional constraints. I think using such a custom distribution / extension should be working without any additional efforts on the Turing side.

[richardreeve]

Okay. The reason I assumed that might not be the case is that the documentation says that to use a custom distribution you need to define logpdf() for that distribution, and I wasn't confident how to do that given that the pdf would just be a rather unhelpful Inf at the Dirac delta?

[trappmartin]

True. I think logpdf is unfortunately used for the log mass function in the Distributions package. Or it’s a alias, not sure. Anyways, we assume logpdf = log pdf in the continuous case and logpdf = log pmf in the discrete. I know this is not nice but for now simplifies the interface.

[trappmartin]

We could also dispatch on the ValueSupport and use logpdf or logpmf respectively. However, I’m not sure how to handle mixed support in such a setting. I suppose this is what you are worried about.

[richardreeve]

Yes, that's exactly it. At the moment there is only (log)pdf, which works for both. I'm proposing to replace pdf with pmf for discrete distributions in Distributions.jl but to alias pdf to pmf for backwards compatibility. Then defining both for mixed distributions and dispatching on ValueSupport would work fine for Turing.jl (though I have to admit I still don't really understand how you can handle switching between spike and slab in the likelihood ratio).

[richardreeve]

So, any thoughts about this? I would really like to have a go at updating Distributions.jl's handling of discrete / continuous / discontinuous distributions before it hits 1.0, but there's no point (to me!) in doing that unless I know what would be usable you all... on which note, I would also like to know how you would (like to?) handle a discrete distribution that has support over any countable set (e.g. of strings, symbols) rather than just integers, as there is also a discussion of how to handle them in Distributions.

[trappmartin]

Hi, we will not be able to handle discrete distributions over any countable set as we internally assume that values of any distribution are sub-typed by Real. But I assume we could relax this if necessary without much performance issues. @mohamed82008 is this correct?

Regarding mixing discrete and continuous distributions in the log joint. I don't see why this would be a problem, but I might easily misunderstand the issue you are raising. For discontinuous distributions (assuming IS based inference) this seems unproblematic too. At least in my opinion.

[richardreeve]

Handling non-real values would be great if it's possible, because there is already one PR - JuliaStats/Distributions.jl#916 - to allow non-numeric discrete distributions, and I would be hoping to generalise that.

As far as mixing discrete and continuous distributions, again great if it already works - as I've said before I'm simply a user of this stuff and I was assuming that you'd be calculating some kind of acceptance ratio in some MH-like world, and I wasn't sure what that even means if the priors didn't have equivalent terms on the top and bottom - one is a mass and the other is a density. However, are you suggesting that you wouldn't use that kind of approach at all and just sample repeatedly from the prior(?) for some kind of particle-based IS approach? So that kind of ratio never gets in to the equations. Isn't that incredibly inefficient not making proposals based on the current state and the prior ratio? Sorry for being so ignorant, I'm just many years out of touch with this stuff obviously...

[trappmartin]

Mixing discrete and continuous RVs in particle-based samplers is not an issue.

Thanks for sharing the PR. Having non-numeric discrete distributions would be really nice for many applications. I'm happy to open a PR on the Turing side once your generalisation is in Distributions.jl. Please keep me in the loop.

[richardreeve]

Will do.

[DoktorMike]

Did this ever materialize and left open or was it perhaps downprioritized? Just curious as i think the spike and slab prior is an interesting and useful beast. 😊

[richardreeve]

I was very keen on the JuliaStats/Distributions.jl PRs that might have helped enable it, but after a few months I ran out of time and it seemed that no-one really wanted to merge it anyway, so I gave up (the final version was JuliaStats/Distributions.jl#951 fwiw). I'd love to have it merged, or this worked on independently of that PR...

[yebai]

I feel this is better done in another package, e.g. Distributions.

[richardreeve]

@yebai - in the unlikely event that I have time to take this back up again in Distributions.jl, it would be good to know what actually needs to be implemented for this to work inside Turing.jl... is there a specification anywhere? I know I knew this a couple of years ago, but I'd just like to get myself back up to date...

[yebai]

in the unlikely event that I have time to take this back up again in Distributions.jl, it would be good to know what actually needs to be implemented for this to work inside Turing.jl

thanks, @richardreeve for the message. Turing is compatible with Distributions API, so there is nothing to do from the Turing side!

[devmotion]

allow discrete distributions have support over any countable set rather than just integers

This is supported by Distributions nowadays, e.g. DiscreteNonParametric and Dirac allow non-integer support. It seems Discrete and DiscreteDistribution are (more or less commonly) understood as a countable support and distributions with countable support (i.e., distributions of discrete random variables in line with e.g. the simple notion on Wikipedia), and Continuous and ContinuousDistribution are viewed to refer to uncountable support and distributions with uncountable support (i.e., distributions of continuous random variables). This seems to deviate from their original interpretation but is consistent with the implemented DiscreteDistributions and ContinuousDistributions, allows to characterize each distribution as either DiscreteDistribution or ContinuousDistribution, and seemed the practically most useful approach for e.g. defining product distributions of discrete and continuous distributions (as e.g. in JuliaStats/Distributions.jl#1391) without having to break or change the type structure.

However, unfortunately multiple (fallback) definitions still incorrectly assume that the support of DiscreteDistributions is integer-valued. In general, I think a better type system for the support of discrete (but also continuous) distributions, possibly partly inspired or taken from e.g. DomainSets, would be very very helpful. This could e.g. also fix some bugs with countably infinite support such as for Poisson.