BayesianTools.jl is a Julia package with methods useful for Monte Carlo Markov Chain simulations. The package has two submodules:

• ProductDistributions: defines a ProductDistribution type and related methods useful for defining and evaluating independent priors
• Link: useful to rescale MC proposals to live in the support of the prior densities

The package is registered

(v1.x) pkg> add BayesianTools

The following code shows how a product distribution resulting from multiplying a normal and a Beta can be obtained

using BayesianTools.ProductDistributions
p = ProductDistribution(Normal(0,1), Beta(1.,1.))
n = length(p) ## 2 -> Number of distributions in the product

To check whether an Array{Float64} is in the support of p

insupport(p, [.1,2.]) ## false
insupport(p, [.1,1.]) ## true

The logpdf and the pdf at a point x::Array{Float64}(n) are

logpdf(p, [.1,.5]) # = logpdf(Normal(0,1), .1) + logpdf(Beta(1.,1.), .5)
pdf(p, [.1,.5]) # = pdf(Normal(0,1), .1) * pdf(Beta(1.,1.), .5)

It is also possible to draw a sample from p

rand!(p, Array{Float64}(2,100))

invlink and link are useful to transform and back-transform the parameters of a parametric statistical model according to the support of its distribution. logjacobian provides the log absolute Jacobian of the inverse transformation applied by invlink.

The typical use case of the methods in the Links is best understood by an example. Suppose interest lies in sampling from a Gamma(2,1) distribution $$\pi(x) = xe^{-x}, \quad x \geq 0 \text{ .}$$

This is a simple distribution and there are many straightforward ways to draw from it. However, we will consider employing a random walk Metropolis-Hastings (MH) sampler with a standard Gaussian proposal.

The support of this distribution is $x>0$ and there are four options regarding the proposal distribution:

1. Use a Normal(0,1) and proceed as you normally would if the support of the density was (-Inf, +Inf)
2. Use a truncated normal distribution
3. Sample from a Normal(0,1) until the draw is positive
4. Re-parametrise the distribution in terms of $y=\exp(y)$ and draw samples from $$\tilde{\pi}(y) = \log(y)e^{-\log(y)} \text{ .}$$

The first strategy will work just fine as long as the density evaluates to 0 for values outside its support. This is the case for the pdf of a Gamma in the Distributions.jl package.

The second and the third strategy are going to work as long as the acceptance ratio includes the normalizing constant (see Darren Wilkinson's post).

The last strategy also needs an adjustment to the acceptance ratio to incorporate the Jacobian of the transformation.

The code below use invlink, link, and logjacobian to carry out the r.v. transformation and the Jacobian adjustment.

Notice that the Improper distribution is a subtype of ContinuousUnivariateDistribution. Links defines methods for Improper that allow the transformations to go through automatically. (Improper can also be used as a component of the ProductDistribution which is useful if an improper prior was elicited for some components of the parameter.)

using BayesianTools.Links
function mcmc_wrong(iters)
   chain = Array{Float64}(iters)
   gamma = Gamma(2, 1)
   d = Improper(0, +Inf)
   lx  = 1.0
   for i in 1:iters
      xs = link(d, lx) + randn()
      lxs = invlink(d, xs)
      a = logpdf(gamma, lxs)-logpdf(gamma, lx)       
      (rand() < exp(a)) && (lx = lxs)
      chain[i] = lx
   end
   return chain
end
function mcmc_right(iters)
   chain = Array{Float64}(iters)
   gamma = Gamma(2, 1)
   d = Improper(0, +Inf)
   lx  = 1.0
   for i in 1:iters
      xs = link(d, lx) + randn()
      lxs = invlink(d, xs)
      a = logpdf(gamma, lxs)-logpdf(gamma, lx)
      ## Log absolute jacobian adjustment
      a = a - logjacobian(d, lxs) + logjacobian(d, lx)
      (rand() < exp(a)) && (lx = lxs)
      chain[i] = lx
   end
   return chain
end

The results are

mc0 = mcmc_wrong(1_000_000)
mc1 = mcmc_right(1_000_000)
using Plots
Plots.histogram([mc0, mc1], normalize=true, bins = 100, fill=:slategray, layout = (1,2), lab = "draws")
title!("Incorrect sampler", subplot = 1)
title!("Correct sampler", subplot = 2)
plot!(x->pdf(Gamma(2,1),x), w = 2.6, color = :darkred, subplot = 1, lab = "Gamma(2,1) density")
plot!(x->pdf(Gamma(2,1),x), w = 2.6, color = :darkred, subplot = 2, lab = "Gamma(2,1) density"))
png("sampler")