PolyaGammaHybridSamplers.jl

PolyaGammaHybridSamplers.jl is a Julia package that implements a hybrid sampling approach for the Pólya-Gamma distribution. The Pólya-Gamma distribution is a continuous distribution on the positive real number line and is most often used for latent variable augmentation in models like logistic regression or logit stick-breaking processes. It is parameterized by two values, $b > 0$ and $z\in\mathbb{R}$, and is often assigned the distributional notation $PG(b,z)$. An excellent tutorial on using the Pólya-Gamma distribution can be found at this blog post by Gregory Gunderson.

Background Details

The class of Pólya-Gamma distributions was first presented in a paper by Polson et al. (2013). There is no closed-form expression for the $PG(b,z)$ density function $f$, but is rather defined through its Laplace transform:

$$ \mathcal{L}_f(t) = \frac{\cosh^b(\frac{z}{2})}{\cosh^b\left(\sqrt{\frac{z^2/2 + t}{2}}\right)} $$

This can be derived from an alternative definition where the $PG(b,z)$ distribution is represented as an infinite mixture of gamma variables. The authors comment that a naive approach to sampling could be to generate atruncated number of gamma values and combine them using the mixture formula, but they warn that this approximation is both inefficient and potentially dangerous.

In order to develop an effecient yet exact sampling method, the authors drew connections to a subclass of $J^*$ distributions discussed in Devroye (2009). The ultimate conclusion is an extremely efficient rejection sampler, called the Devroye method, that is used to simulate from a $PG(1,z)$ distribution. If $b$ is an integer, then the $PG(b,z)$ distribution can be drawn by summing $b$ independent draws from the $PG(1,z)$ distribution.

The following year, a technical report written by Windle et al. (2014) was published on arXiv detailing an alternative sampling technique and an approximate technique. The alternative technique requires $b\in[1,4]$, but allows for otherwise non-integer $b$ values. The approximate algorithm, called the saddle-point sampler, involves approximating the density using the moment generating function and then constructing a smart rejection sampler. The saddle-point sampler is both fast and increasingly accurate for larger values of $b$. The report concludes with a performance comparison between the naive gamma mixture, the Devroye method, the alternative technique, the saddle-point sampler, and a normal approximation. They subsequently suggest the following hybrid sampler scheme to balance accuracy and efficiency:

Sampler	Ideal Circumstances
Devroye Method	$b = 1,2$
Alternative Technique	$3 \le b \le 13$ and non-integer $1 < b < 3$
Saddle-point Sampler	$13 < b < 170$
Normal Approximation	$b \ge 170$

Previous Code Implementations

The authors of Polson et al. (2013) published the BayesLogit R package with C++ code implementing the Devroye sampling method, later expanding it to include the techniques presented in Windle et al. (2014). The source code is available at the BayesLogit GitHub site. However, it is important to note that, as of version 2.1, inline comments for the main sampling function, rpg_hybrid, state that the alternative technique code needs further review and/or development (source, lines 138-141) and instead the naive gamma mixture is used for small, non-integer $b$ values combined with the Devroye method for integer $b>1$.

Interest in the Pólya-Gamma distribution within the Julia community appears first in a GitHub issues thread started by a user called currymj under the main Distributions.jl repository. Ultimately, currymj developed a separate package called PolyaGammaDistribution.jl that was essentially a direct translation of the BayesLogit Devroye method code for integer $b$ parameters.

However, Julia would eventually update to v1.0 and begin requiring a different package structure. When asked about future compatibility, a collaborator, user maximerischard, indicated in a GitHub issues thread that there were no plans to upgrade the package. The user, igutierrezm, who started the thread would then go on to develop a Julia v1.0 compatible package, PolyaGammaSamplers.jl. This new package implemented the Devroye method similar to its predecessor, but using the new sampler framework, which is a more limited implementation compared to the full distribution requirements of Distributions.jl, but is also more fitting for the Pólya-Gamma distribution.

Over time, a handful of variations came about, including a PolyaGammaDistribution.jl fork by user yunzli (link) and a rewritten Devroye method sampler in AugmentedGaussianProcesses.jl by user theogf. Perhaps the most recent discourse comes from another GitHub issues thread, where user theogf asked whether their Devroye method code should be merged into Distributions.jl. The issue was ultimately never resolved due to conflicts between the MIT and GPL-3.0 licenses.

Scope of PolyaGammaHybridSamplers.jl

As the name suggests, the PolyaGammaHybridSamplers.jl package attempts to fill a hole present in previously developed Julia packages by implementing other sampling techniques for the Pólya-Gamma distribution beyond the Devroye method. It centers around a hybrid sampler which selects the most appropriate sampling technique based on the value of $b$. The hybrid sampler is similar to the one given by Windle et al. (2014):

Sampler	Condition
Normal Approximation	$b \ge 170$
Saddle-point Sampler	$13 \le b < 170$
Devroye Method	$b = 1,2,...,13$
Gamma Mixture	$0 < 1 < b$
Devroye + Gamma Mixture	$1 < b < 13$ non-integer $b$

To generate a random variate, create a sampler object and then use rand (vectorized versions of rand and rand! are also available):

pg = PolyaGammaHybridSampler(4, 1.5)
data = rand(pg)

The type of sampler can controlled by including, for example, SADDLEPOINT as an argument:

pg = PolyaGammaHybridSampler(4, 1.5, SADDLEPOINT)
data = rand(pg)

Other valid options include DEVROYE, NORMALAPPROX, GAMMAMIXTURE, and DEVROYEGAMMAMIXTURE. But be careful as no warning about the approximation quality or efficiency will be given if the sampler is forced to use a technique that is not appropriate for the given parameters.

License Issue and Merging with Distributions.jl

This package is largly a translation of the original Polson et al. (2013) C++ code in the package BayesLogit, including the later developments described in Windle et al. (2014). Some structural inspiration is also taken from the Julia packages PolyaGammaDistribution.jl (currymj), PolyaGammaSamplers.jl (igutierrezm), and AugmentedGaussianProcesses.jl (theogf). All of these packages are licensed under the GPL-3.0 license, except for AugmentedGaussianProcesses.jl which is licensed under the MIT "Expat" license. Therefore, the PolyaGammaHybridSamplers.jl package is necessarily licensed under GPL-3.0. The question of whether this code can merged into Distributions.jl has likely been answered in this GitHub issues thread started by theogf: the MIT license of Julia's base environment conflicts with the GPL-3.0 requirements, so probably not.