BinomialLikelihood for count data / bucketed BernoulliLikelihoods

[feynmanliang]

Hi, I have two questions related to a potential application of gpytorch for Binomial regression of count data.

1. Does there exist a way to specify a Binomial likelihood (eg p(trials, pos | f) GP(f) = Binomial(pos | prob=Phi(f), trials) GP(f)) and use gpytorch to Monte-Carlo approximate the marginal likelihood / perform Bayesian inference over GP hyperparameters?
2. Is there a reason why ExactGP is restricted to only Gaussian likelihoods? I understand we may no longer be able to do Schur complement tricks, but if one is interested in Bayesian inference of hyperparameters (eg https://mc-stan.org/docs/stan-users-guide/fit-gp.html#logistic-gaussian-process-regression) then is it still possible to exactly / closely approximate kernels terms required for computing the joint log probability?

[gpleiss]

1. Does there exist a way to specify a Binomial likelihood (eg p(trials, pos | f) GP(f) = Binomial(pos | prob=Phi(f), trials) GP(f)) and use gpytorch to Monte-Carlo approximate the marginal likelihood / perform Bayesian inference over GP hyperparameters?

This would require a custom likelihood, but I imagine it shouldn't be too challenging to write one. See the BetaLikelihood as an example.

(This is something that would be valuable in GPyTorch, so if you do write one, then please submit a PR. It'd also be great to get feedback/doc suggestions for implementing a custom likelihood.)

[gpleiss]

@danbider I know you were looking for a good first PR - if you're still interested this would be a good issue!

[feynmanliang]

ATM I am able to achieve the desired effect by delegating the likelihood the the enclosing PPL. e.g. to extend https://docs.gpytorch.ai/en/stable/examples/01_Exact_GPs/GP_Regression_Fully_Bayesian.html to Binomial data (x=covariate, n=total trials, y=num successes):

def pyro_model(x, y):
    with gpytorch.settings.fast_computations(False, False, False):
        sampled_model = gp_model.pyro_sample_from_prior()
        f_x = pyro.sample("f_x", sampled_model.likelihood(sampled_model(x)))
        pyro.sample("obs", Binomial(total_count=trials, probs=dist.Normal(0, 1).cdf(f_x)), obs=y)
    return y

Effectively the sample performs a Monte-Carlo approximation to f_x. I think this is appropriate: if an end user wants an explicit Monte-Carlo sampling approximation, then they should use the sample API in their implementing PPL. Thoughts?

IIUC this is different than a BinomialLikelihood, which we would behave similarly to BetaLikelihood and used with a VariationalGP to define a variational approximation objective?

[gpleiss]

That's sort of what I'm thinking for an exact GP MCMC interface. But here's a few modifications to your existing code that I think should be necessary for it to work:

gp_model.eval()  # so you don't get the `you must train on training inputs` error message.

def pyro_model(x, y):
    with gpytorch.settings.fast_computations(False, False, False):
        sampled_model = gp_model.pyro_sample_from_prior()
        with gpytorch.settings.prior_mode():  # ensure that sampled_model(x) returns the prior dist
            f_x = pyro.sample("f_x", sampled_model(x))  # no likelihood - that'll add i.i.d. gaussian noise
        pyro.sample("obs", Binomial(total_count=trials, probs=dist.Normal(0, 1).cdf(f_x)), obs=y)
    return y

IIUC this is different than a BinomialLikelihood, which we would behave similarly to BetaLikelihood and used with a VariationalGP to define a variational approximation objective?

pyro.sample("obs", Binomial(total_count=trials, probs=dist.Normal(0, 1).cdf(f_x)), obs=y) is essentially what a BinomialLikelihood is doing. I imagine that we could have two Pyro interfaces (similar to what we have for variational models):

1. (High level) If you are trying to build a standard predictive model, then you can define a GPyTorch GP + likelihood, and we handle the pyro stuff under the hood.
2. (Low level) alternatively, you could define a pyro model like the one you have: where you are actually making the pyro.sample call yourself for the observational model.

[gpleiss]

BinomialLikelihood PR discussions ⬇️

[danbider]

a strategy (VI?) for posterior inference and hyperparam learning. But I now saw you attached pyro's variational models API?

[gpleiss]

Follow this PR for guidelines: #1211

The one thing that this PR missed was adding the new likelihoods to the documentation, and adding docstrings for all the parameters. See:

• 990640f#diff-ebcbdab85136c737a7c228222b840bf5c5da5a030946a7a13d5fd56fbafcb02d
• 990640f#diff-b9dfd3f0f8bf0360ba9765110aaa62ba185b35cc7da4903080207f1cbabbad06

It should be as simple as:

• Adding a BinomialLikelihood module
• Adding a unit test file for BinomialLikelihood
• Adding BinomialLikelihood to the docs

[gpleiss]

a strategy (VI?) for posterior inference and hyperparam learning. But I now saw you attached pyro's variational models API?

That already exists. The variational inference modules are likelihood agnostic.

[Balandat]

@mshvartsman, @crasanders did you implement this at some point?

[mshvartsman]

We haven't. We've got Bernoulli with non-Probit links (https://github.com/facebookresearch/aepsych/blob/main/aepsych/likelihoods.py#L14, links at https://github.com/facebookresearch/aepsych/blob/main/aepsych/acquisition/objective.py) and soon a categorical likelihood for ordinal models (facebookresearch/aepsych#152), and the MC expectations all seem to happen under the hood just fine (at least for the ELBO, we haven't tried any other MLLs).

We have had a bunch of shape errors trying to use the ordinal likelihood with tests subclassing from the gpytorch likelihood test base classes though.

[gpleiss]

2. Is there a reason why ExactGP is restricted to only Gaussian likelihoods? I understand we may no longer be able to do Schur complement tricks, but if one is interested in Bayesian inference of hyperparameters (eg https://mc-stan.org/docs/stan-users-guide/fit-gp.html#logistic-gaussian-process-regression) then is it still possible to exactly / closely approximate kernels terms required for computing the joint log probability?

This is something that we're very much interested in implementing. I think the API design might take a bit of thought, but I think it could be something similar to the low level Pyro API for variational models.

[gpleiss]

(Moved this discussion to a new thread, since the other thread accidentally got overtaken with BinomialLikelihood PR discussions 😅)