add helper functions for marginalising discrete random variables

[goldingn]

greta can't currently sample discrete random variables, because HMC is the only implemented sampler (though see #104)

You can use HMC to sample (some) models with latent discrete random variables, by explicitly marginalising out the random variable. This is possible in the dev version of greta, using the new mixture() distribution.

It would be great to provide utility functions to help people do that marginalisation though. The computation for these models isn't too difficult, but finding a nice user interface could be a bit tricky.

This gist: https://gist.github.com/goldingn/4192ddff71037ca3f18927b8d273cd2a has a prototype implementation for a latent Poisson random variable in a simple model. That interface requires the user to pass a function for the model likelihood conditional on the scalar Poisson variable N, a greta array lambda for the rate parameter, and the maximum value of N up to which to integrate the likelihood.

The syntax looks like this:

likelihood <- function (N) bernoulli(0.9 ^ N)
distribution(y) <- marginal_poisson(likelihood,
                                    lambda,
                                    max_n = 10)

The gist shows how to get the posterior over the weights for different values of N, and how to sample values of N post-hoc.

A more general interface would be nice, i.e. just marginalise(), but with a discrete greta distribution over which to marginalise, e.g.:

distribution(y) <- marginalise(likelihood,
                               poisson(lambda),
                               values = 1:10)

[goldingn]

The more general approach could use the distribution node to directly calculate the log PMF at the values, which would be more stable than this prototype (if they can be passed to mixture())

[goldingn]

The pattern for this should be:

marginalise(fun,
            distribution,
            ...,
            method = discrete())

and used like this:

lambda <- lognormal(0, 1)
theta <- normal(0, 1)
sd <- lognormal(0, 1)
y <- rnorm(100)

# # with a sampleable discrete RV the rest would be:
# n <- poisson(lambda)
# mu <- theta ^ n
# distribution(y) <- normal(mu, sd)

# instead:
fun <- function (n, theta, sd) {
  mu <- theta ^ n
  distribution(y) <- normal(mu, sd)
}

marginalise(fun,
            poisson(lambda),
            theta = theta,
            sd = sd,
            method = discrete(values = 0:10))

This interface enable extension to other forms of marginalisation (e.g. over continuous distributions) which is planned for a future release

[goldingn]

There's a prototype marginalisation interface in the feature/marginalisation branch. It handles both discrete variables, and variables with multivariate normal priors (with the Laplace approximation). I'd like to add more rigorous tests of the Laplace approximation before I merge.

I also need to edit this to return relevant summary statistics of the marginalised variables; posterior samples for the weights (in the discrete case) or the mode and covariance (in the Laplace approximation case). This should be feasible, by returning that as the output of the marginalisation

The interface is as described above, however it should be possible to implement the following interface, which is much nicer, and would lead to much more useful model plots:

# marginalise a discrete random variable and carry out HMC on the
# continouous variables:
y <- rnorm(100)

# priors
lambda <- lognormal(0, 1)
theta <- normal(0, 1)
sd <- lognormal(0, 1)

# model
n <- poisson(lambda)
mu <- theta ^ n

# likelihood
distribution(y) <- normal(mu, sd)

# n is a discrete random variable so can't be sampled by HMC, but we can
# marginalise it with the following line:
marginalise(n, method = discrete_marginalisation(0:10))

m <- model(lambda)
draws <- mcmc(m, hmc())

to implement this interface, we would need to:

• flag n for marginalisation in marginalise() (add a flag in the node)
• construct the dag as normal (use this for plotting, with the marginalised variable indicated)
• find the nodes that depend on n and turn them into a conditional density function (this is the technical bit)
• find the distribution for n, and its parameters
• check that function is valid, as before
• construct a second, marginalised dag and use that for computation

all of that is very doable, if somewhat time consuming

[jeffreypullin]

Hi Nick,

I'm (coincidentally!) currently doing a project on bayesian annotation models with categorical predictors (lots of discrete parameters) - at this stage just getting my head round how to write them in Stan so I'd love to help/ contribute some example models when this is done.