Feature request: Distributions for Jeffreys-type priors

It would be great if nimble could use Jeffreys-type improper distributions for normal and multivariate normal models. Nimble could exploit the conjugacy of these priors.

For example, if $y \sim N(\mu, \sigma^2)$ , the Jeffreys prior is $\pi(\mu, \sigma^2) \propto \left( \sigma^2 \right)^{-1/2}$ . More generally, it would be useful to have a prior of the form $\pi(\mu, \sigma^2 | d) \propto \left( \sigma^2 \right)^{-d/2}$ . For example, in a linear regression setting, fixing d = p (number of regression coefficients) yields similar results as a frequentist analysis. Setting d = 0 yields a uniform improper prior (the flat prior), so that this is a flexible generalization.

For the multivariate normal model, it would be useful to have priors of the form $\pi(\mu, \Sigma | d) \propto \left\lvert \Sigma \right\rvert^{-d/2}$ . For example, d= p+2 corresponds to Jeffreys prior, d = p corresponds to objective Bayes. Of course, d = 0 corresponds to a uniform prior across the space of positive definite covariance matrices. Thus, these priors are similar to LKJ priors for correlation matrices, except that they are for covariance matrices and are conjugate for normal models.

Although one could easily code these priors in Nimble, I do not believe one can make Nimble recognize the conjugacy. Here's my attempt at creating the prior for the multivariate normal case

djeffprior <- nimbleFunction(
  run = function(
    x = double(2)
    , pow = double(0, default = 1)
    , log = integer(0, default = 0)
  ) {
    returnType(double(0))
    logProb <- -1.0 * pow * logdet(x)
    if(log) return(logProb)
    else return(exp(logProb))
  }
)

Configuring the MCMC, it uses a RW_block sampler for Sigma:

===== Monitors =====
thin = 1: gamma, mu, Sigma
===== Samplers =====
RW_block sampler (1)
  - Sigma[1:2, 1:2] 
conjugate sampler (3)
  - gamma
  - mu[1:2, 1]
  - mu[1:2, 2] 
binary sampler (100)
  - c0[]  (100 elements)

Hi @ethan-alt

I don't have insights on making this a possible feature (I'll leave the developers to chime in on this), however I can share a possible workaround in case it is useful: a "trick" used in old BUGS examples for linear regression is to approximate the Jefferey's prior for the inverse of the variance using a Gamma distribution with very small values for the parameters (e.g. tau ~ dgamma(0.001, 0.001), where tau = 1/sigma^2).

I don't have a better reference than these these slides at the moment (see page 122). I did look into this some time ago for linear regression, but there might be similar workarounds for other cases.

Alternatively, one can write a custom sampler for Sigma, since the posterior will be known and proper. In the user manual there is an example of a custom sampler in case is something you may be interested in (see https://r-nimble.org/html_manual/cha-progr-with-models.html#sec:user-samplers ).

@ethan-alt Great question, and this is one aspect of NIMBLE that I actually have some knowledge of.

As you noted, writing and using a custom-writtten prior distribution for the covariance matrix Sigma is straightforward.

Having nimble automatically detect and implement a conjugate sampler for it is a little more effort, but not impossible. Actually, it would be quite straightforward to put into the code base of nimble, and it would work "out of the box" for you. What would be more difficult is for you to load nimble (with the currently implemented conjugacy relationships), then modifying the nimble internals (specifically, of the conjugacy system) to inform it about the new conjugacy, and have it automatically detect it. Nimble was written to very easily add new (custom written) distributions, but the conjugacy system was not written with the idea of users defining and adding their own (new) conjugate relationships between distributions.

If you want more insight about how conjugacies are defined inside nimble, I suggest you look at: https://github.com/nimble-dev/nimble/blob/devel/packages/nimble/R/MCMC_conjugacy.R

specifically the list object defined at the very top, conjugacyRelationshipsInputList, which is the complete definition of all conjugate relationships used by nimble. The conjugate samplers themselves are dynamically generated (as-needed, on the fly) from the information in this list. (So, adding this request to nimble would require (1) adding the prior distribution itself, then (2) adding an entry to this list to define the conjugate relationship).

I would invite @paciorek thought's on this prior distributions and the conjugacies you've raised, and whether we might consider adding them to the package. (Noting also that we do not currently have a vectorized version of the dflat distribution, which I believe is what you'd want for the mean vector "mu" for the MVN conjugacy.)

@ethan-alt To get you going, the most direct path would be :

Writing your prior distribution, as a nimbleFunction, as you already have.
Writing the relevant conjugate sampler yourself (we can provide guidance if that's helpful).
In your R code, creating the default MCMC, then removing nimble's default sampler, and adding your custom-written conjugate sampler, which is extremely straight-forward by modifying the MCMC configuration object.

Let me know if this helps. And great question.

For the univariate normal case, you should be able to get what you want with:

y ~ dnorm(mu, var = sigma2)
mu  ~ dflat()
sigma2 ~ invgamma(-0.5, 0)

You should get conjugate samplers assigned by default.

For the MVN case, we'd need to do a bit of work. Not much, but lots of competing priorities for development time...

nimble-dev / nimble

Feature request: Distributions for Jeffreys-type priors #1196