## BayesianTools.jl

Julia package designed to handle product distributions
Author gragusa
Popularity
11 Stars
Updated Last
1 Year Ago
Started In
September 2016

# BayesianTools.jl

`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

## Installation

The package is registered

`(v1.x) pkg> add BayesianTools`

## Usage

### ProductDistributions

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

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 and draw samples from

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` package.

The second and the third strategy is 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()
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()
a = logpdf(gamma, lxs)-logpdf(gamma, lx)
a = a - logjacobian(d, lxs) + logjacobian(d, lx)
(rand() < exp(a)) && (lx = lxs)
chain[i] = lx
end
return chain
end```

The results is

```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")```

### Used By Packages

No packages found.