Gamma prior for precision parameter

[osorensen]

In the JMLR paper, an exponential prior was assumed for the precision parameter. This makes it hard to obtain a "weakly informative" prior, as the only option we have is to concentrate mass close to 0. In contrast, Marta's 2019 AAS paper used a gamma prior, which is much more flexible. Here are some examples, the first two plots showing the exponential distribution and the second two plots showing the gamma distribution. If we have reason to believe the alpha is away from zero but also much smaller to 1000, then the last plot with the gamma distribution is definitely the best thing to use.

x <- seq(from = 0, to = 40, by = .1)

plot_expdist <- function(lambda) {
  par(mfrow = c(2, 2))
  plot(x, pexp(x, rate = lambda), xlab = expression(alpha), type = "l")
  plot(x, dexp(x, rate = lambda), xlab = expression(alpha), type = "l")
  s <- rexp(1000, rate = lambda)
  hist(s, xlab = expression(alpha))
  plot(density(s), xlab = expression(alpha))
}

plot_expdist(.001)

plot_expdist(.3)

plot_gammadist <- function(gamma, lambda) {
  par(mfrow = c(2, 2))
  plot(x, pgamma(x, shape = gamma, rate = lambda), xlab = expression(alpha), type = "l")
  plot(x, dgamma(x, shape = gamma, rate = lambda), xlab = expression(alpha), type = "l")
  s <- rgamma(1000, shape = gamma, rate = lambda)
  hist(s, xlab = expression(alpha))
  plot(density(s), xlab = expression(alpha))
}

plot_gammadist(1, .001)

plot_gammadist(2, .1)

^{Created on 2024-02-07 with reprex v2.1.0}

Also, this does not need to create any breaking changes, since the exponential distribution is a special case of the gamma distribution, as shown in the example below:

x <- seq(from = 0, to = 40, by = .1)
plot(x, dexp(x), type = "o")
lines(x, dgamma(x, 1), col = 2)

^{Created on 2024-02-07 with reprex v2.1.0}

[osorensen]

This means that in the ratio here:

https://github.com/ocbe-uio/BayesMallows/blob/a1d3ee857fdfb5fc0a3934c82db5507788fbeecb/src/proposal_functions.cpp#L21-L28

I need to add a term (gamma - 1) * std::log(alpha_proposal / alpha_old).

[osorensen]

Actually, two terms cancel and it becomes

 double ratio = 
   alpha_diff / n_items * rank_dist + 
   priors.lambda * alpha_diff + 
   sum(observation_frequency) * 
   (pfun->logz(alpha_old) - pfun->logz(alpha_proposal)) + 
   priors.gamma * (std::log(alpha_proposal) - std::log(alpha_old));