Data Augementation with Polya Inverse Gamma

TL;DR

This paper introduces a novel data augmentation method using Pólya Inverse Gamma distributions to improve inference in models involving gamma functions, applicable to various statistical and machine learning models.

Contribution

It develops the theory of Pólya Inverse Gamma distributions and demonstrates scalable EM and MCMC algorithms for models with gamma functions, addressing a key inference challenge.

Findings

Effective data augmentation for gamma-related models

Scalable algorithms for inference in complex models

Successful application to negative binomial regression and Dirichlet allocation

Abstract

We use the theory of normal variance-mean mixtures to derive a data augmentation scheme for models that include gamma functions. Our methodology applies to many situations in statistics and machine learning, including Multinomial-Dirichlet distributions, Negative binomial regression, Poisson-Gamma hierarchical models, Extreme value models, to name but a few. All of those models include a gamma function which does not admit a natural conjugate prior distribution providing a significant challenge to inference and prediction. To provide a data augmentation strategy, we construct and develop the theory of the class of P\'olya Inverse Gamma distributions. This allows scalable EM and MCMC algorithms to be developed. We illustrate our methodology on a number of examples, including gamma shape inference, negative binomial regression and Dirichlet allocation. Finally, we conclude with directions for future research.

Tables

Table 1. Table 1: List of Gamma Models

Distribution	Likelihood	Applications
Gamma	$\frac{{\beta }^{\alpha }}{\Gamma \left(\alpha \right)}{x}^{\alpha -1}{e}^{-\beta x}$	Gamma process, Poisson regression
Inverse Gamma	$\frac{{\beta }^{\alpha }}{\Gamma \left(\alpha \right)}{x}^{-\alpha -1}{e}^{-\beta /x}$	survival analysis, conjugate prior
Beta	$\frac{\Gamma \left(\alpha +\beta \right)}{\Gamma \left(\alpha \right)\Gamma \left(\beta \right)}{x}^{\alpha -1}{\left(1-x\right)}^{\beta -1}$	order statistics, wavelet analysis
Dirichlet	$\frac{\Gamma \left({\sum }_{k=1}^K{\alpha }_k\right)}{{\prod }_{k=1}^K\Gamma \left({\alpha }_k\right)}{\prod }_{k=1}^K{x}_k^{{\alpha }_k-1}$	topic model, Bayesian networks
Negative Binomial	$\frac{\Gamma \left(x+\alpha \right)}{\Gamma \left(\alpha \right)\Gamma \left(\alpha +1\right)}{p}^r{\left(1-p\right)}^x$	stock control problems, negative binomial regression

Table 2. Table 2: Integral Representation for Gamma Functions

	Integral Representation	Auxiliary Variables
$\Gamma \left(\alpha \right)$	${\int }_0^{\infty }{x}^{\alpha -1}{e}^{-x}𝑑x$	Gamma
$\frac{1}{\Gamma \left(\alpha \right)}$	${\int }_0^{\infty }\alpha e^{-{\alpha }^{2}x+\gamma \alpha }{p}_0\left(x\right)𝑑x$	P-IG
$\frac{\Gamma \left(\alpha \right)}{\Gamma \left(\alpha +\beta \right)}$	${\int }_0^1\frac{1}{\Gamma \left(\beta \right)}{x}^{\alpha -1}{\left(1-x\right)}^{\beta -1}𝑑x$	Beta

Table 3. Table 3: Posterior Moments

$n$	$x_a$	$x_g$	Method	Mean	Variance	Skewness	Kurtosis
5	7.19	6.05	Numerical Integration	4.768	5.399	0.997	4.494
			Posterior Sampling	4.779	5.630	1.305	6.148
10	5.57	5.01	Numerical Integration	6.271	5.780	0.783	3.921
			Posterior Sampling	6.101	4.857	0.634	3.427
30	5.09	4.26	Numerical Integration	3.252	0.585	0.490	3.361
			Posterior Sampling	3.250	0.605	0.537	3.155

Table 4. Table 4: Homogeneous Setting

		$b_0=0.1$				$b_0=1$				$b_0=5$
		Root Mean Square Error ($\times {10}^3$)
$S$	$K$	DA	DA-N	DA-E	MH	DA	DA-N	DA-E	MH	DA	DA-N	DA-E	MH
100	10	20.93	20.95	18.95	20.56	20.25	20.36	18.31	19.99	20.22	20.25	18.38	20.19
100	50	7.78	8.03	6.93	5.88	7.71	7.95	6.86	6.03	7.66	7.95	6.84	5.95
1000	10	6.54	6.53	5.95	7.13	6.58	6.54	6.01	7.26	6.38	6.36	5.80	6.87
1000	50	2.49	2.51	2.19	2.48	2.49	2.51	2.20	2.50	2.48	2.50	2.18	2.46
		Effective Sample Size Ratio
$S$	$K$	DA	DA-N	DA-E	MH	DA	DA-N	DA-E	MH	DA	DA-N	DA-E	MH
100	10	0.32	0.33	0.32	0.04	0.32	0.32	0.32	0.04	0.32	0.32	0.32	0.04
100	50	0.08	0.08	0.07	0.01	0.08	0.08	0.07	0.01	0.08	0.08	0.07	0.01
1000	10	0.32	0.32	0.32	0.01	0.33	0.33	0.32	0.01	0.32	0.32	0.32	0.01
1000	50	0.08	0.08	0.07	0.01	0.08	0.07	0.07	0.01	0.08	0.08	0.07	0.01
		Running Time
$S$	$K$	DA	DA-N	DA-E	MH	DA	DA-N	DA-E	MH	DA	DA-N	DA-E	MH
100	10	32.76	0.63	0.22	0.93	31.52	0.61	0.22	0.89	33.37	0.62	0.23	0.89
100	50	174.29	3.61	1.12	2.51	166.49	3.43	1.08	2.43	165.67	3.43	1.07	2.45
1000	10	236.33	5.13	1.31	9.47	206.65	4.59	1.21	8.70	204.72	4.62	1.22	8.79
1000	50	1099.74	30.67	6.18	26.05	1012.35	28.12	5.77	23.99	1014.13	28.14	5.75	23.90

Table 5. Table 5: Heterogeneous Setting

		$b_0=0.1$				$b_0=1$				$b_0=5$
		Root Mean Square Error ($\times {10}^3$)
$S$	$K$	DA	DA-N	DA-E	MH	DA	DA-N	DA-E	MH	DA	DA-N	DA-E	MH
100	10	82.76	82.25	76.54	82.34	82.56	81.90	77.33	83.49	81.95	82.11	78.46	86.96
100	50	81.50	81.41	75.83	553.05	78.53	78.39	73.47	553.83	81.99	81.98	79.15	553.38
1000	10	26.11	26.00	24.50	41.00	25.76	25.61	24.18	41.08	25.88	25.71	24.25	41.39
1000	50	24.80	24.79	23.19	553.35	25.22	25.11	23.66	553.30	25.26	25.21	23.71	552.86
		Effective Sample Size Ratio
$S$	$K$	DA	DA-N	DA-E	MH	DA	DA-N	DA-E	MH	DA	DA-N	DA-E	MH
100	10	0.47	0.49	0.50	0.03	0.48	0.51	0.49	0.03	0.49	0.51	0.51	0.03
100	50	0.44	0.45	0.44	0.01	0.44	0.45	0.44	0.01	0.45	0.46	0.46	0.01
1000	10	0.48	0.50	0.50	0.02	0.47	0.50	0.50	0.01	0.47	0.50	0.50	0.01
1000	50	0.44	0.44	0.44	0.01	0.44	0.45	0.45	0.01	0.44	0.44	0.45	0.01
		Running Time
$S$	$K$	DA	DA-N	DA-E	MH	DA	DA-N	DA-E	MH	DA	DA-N	DA-E	MH
100	10	42.28	0.58	0.26	0.93	32.97	0.49	0.23	0.85	33.36	0.49	0.24	0.86
100	50	186.00	3.07	1.23	2.48	164.91	2.75	1.11	2.25	164.60	2.75	1.11	2.26
1000	10	230.54	3.69	1.27	8.92	205.64	3.33	1.18	8.41	204.56	3.34	1.18	8.51
1000	50	1104.85	22.43	6.25	24.09	1014.25	21.06	5.85	22.41	1013.74	21.06	5.86	22.46

Equations

$$E\left(e^{-s^2{X}_0}\right)={\int }_0^{\infty }{e}^{-s^{2}x}{p}_0(x)dx=\frac{e^{-\gamma s}}{\Gamma (1+s)},s>0.$$

$$\frac{e^{-\gamma s}}{\Gamma (1+s)}=\prod _{k=1}^{\infty }\left(1+\frac{s}{k}\right)e^{-\frac{s}{k}},s\in \mathbb{C}/\{0,-1,-2,...\}$$

$$E\left(e^{-\frac{1}{2}{s}^2{H}_a^{\Gamma }}\right)={\int }_0^{\infty }{e}^{-\frac{1}{2}{s}^{2}x}{p}_a^{\Gamma }(x)dx=\frac{\Gamma (a)}{\Gamma (a+s)}{e}^{\psi (a)s}$$

$$\frac{\Gamma (a)}{\Gamma \left(a+s\right)}{e}^{\psi (a)s}=\prod _{k=0}^{\infty }\left(1+\frac{s}{a+k}\right)e^{-\frac{s}{a+k}}.$$

$$\text{P-IG}(0)\stackrel{D}{=}\frac{1}{2}{H}_1^{\Gamma }.$$

$$E\left(e^{-s^2{X}_{n,0}}\right)={\int }_0^{\infty }{e}^{-s^{2}x}{p}_{n,0}(x)dx=\frac{e^{-\gamma ns}}{\Gamma (1+s{)}^n},n>0,s>0.$$

$$p_c\left(x\right)=Z_c\cdot \exp \left(-c^{2}x\right)p_0(x),x,c>0.$$

$$E\left(e^{-s^2{X}_c}\right)$$

$$=\frac{\Gamma (1+c)}{\Gamma (1+\sqrt{s^2+c^2})}{e}^{-\gamma (\sqrt{s^2+c^2}-c)}.$$

$$\text{P-IG}(c)\stackrel{D}{=}\sum _{k=1}^{\infty }\text{GIG}\left(-\frac{3}{2},2c^2,\frac{1}{2k^2}\right).$$

$$\text{P-IG}(0)\stackrel{D}{=}\sum _{k=1}^{\infty }\frac{1}{4k^2}{\Gamma }_k^{-1}.$$

$$E\left(\sum _{k=N}^{\infty }{G}_k\right)$$

$$\text{Var}\left(\sum _{k=N}^{\infty }{G}_k\right)$$

$$E(g_k)=\frac{1}{2}\left(\frac{1}{k^2+ck}\right),\text{Var}(g_k)=\frac{1}{4c}\left(\frac{1}{k^3+2c{k}^2+c^{2}k}\right).$$

$$E(e^{-s{X}_c})=\exp \left\{{\int }_0^{\infty }(e^{-sx}-1)\nu (x)dx\right\},s\ge 0.$$

$$\nu (x)$$

$$\mu (t)$$

$$\frac{e^{ax}}{\Gamma (1+x)}={\int }_0^{\infty }\varphi (x\mid \mu (\omega ),{\sigma }^2(\omega ))\cdot \tau (\omega )p_0(\omega )d\omega$$

$$p(x)={\int }_0^{\infty }f(x\mid \theta (\omega ))p(\omega )d\omega .$$

$$\text{PG}(b,c)\stackrel{D}{=}\frac{1}{2{\pi }^2}\sum _{k=1}^{\infty }\frac{\text{GIG}(b,2,0)}{(k-1/2{)}^2+c^2/(4{\pi }^2)}.$$

$$\frac{(e^z{)}^a}{(1+e^z{)}^b}={\int }_0^{\infty }f(z\mid \omega )p_{\text{PG}}(\omega \mid b,0)d\omega .$$

$$p(x)=C\cdot \left(\prod _{\ell =1}^L\Gamma (g_{\ell }(x))\right)\left(\prod _{m=1}^M\frac{1}{\Gamma (h_m(x))}\right)\left(\prod _{n=1}^N\frac{\Gamma (j_n(x))}{\Gamma (j_n(x)+{\beta }_n)}\right)\cdot x^{p-1}{e}^{-a{x}^2+bx},x>0,$$

$$p(x)={\int }_{\Omega }p(x,\mathbf{\omega })d\mathbf{\omega }$$

$$p(x)=C\cdot {\int }_{(0,\infty {)}^{L+M}\times (0,1{)}^N}G(x,\mathbf{\tau })\cdot H(x,\mathbf{\omega })\cdot J(x,\mathbf{\eta })\cdot Q(x)\cdot e^{-a{x}^2+bx}d\mathbf{\tau }d\mathbf{\omega }d\mathbf{\eta }$$

$$G(x,\mathbf{\tau })$$

$$H(x,\mathbf{\omega })$$

$$J(x,\mathbf{\eta })$$

$$Q(x)$$

$$p({\tau }_{\ell }\mid x,{\mathbf{\tau }}_{-\ell },\mathbf{\omega },\mathbf{\eta })={\tau }_{\ell }^{g_{\ell }(x)-1}{e}^{-{\tau }_{\ell }}\cdot \frac{\left({\prod }_{{\ell }^{\prime }\ne \ell }{\tau }_{{\ell }^{\prime }}^{g_{{\ell }^{\prime }}(x)-1}{e}^{-{\tau }_{\ell }^{\prime }}\cdot H(x,\mathbf{\omega })\cdot J(x,\mathbf{\eta })\cdot Q(x)\right)}{\int G(x,\mathbf{\tau })\cdot H(x,\mathbf{\omega })\cdot J(x,\mathbf{\eta })\cdot Q(x)d{\tau }_{\ell }},$$

$$p(x)=C\cdot x^{p-1}{e}^{-a{x}^2+bx},x>0$$

$${\tau }_{\ell }\mid x,{\mathbf{\tau }}_{-\ell },\mathbf{\omega },\mathbf{\eta }$$

$${\omega }_m\mid x,\mathbf{\tau },{\mathbf{\omega }}_{-m},\mathbf{\eta }$$

$${\eta }_n\mid x,\mathbf{\tau },\mathbf{\omega },{\mathbf{\eta }}_{-n}$$

$$x\mid \mathbf{\tau },\mathbf{\omega },\mathbf{\eta }$$

$$\tilde{a}$$

Full text

Key Words: Data Augmentation, Pólya Inverse Gamma, Pólya Gamma, Latent Dirichlet Allocation, Gamma Shape, Markov Chain Monte Carlo, Expectation-Maximization;

1 Introduction

Statistical models involving gamma functions are prevalent in statistics and machine learning. For example, topic models, negative binomial regression, time series count models, Poisson-Gamma hierarchical models, non-parametric Bayes, to name but a few (Rossell, 2009; Aktekin, Polson, and Soyer, 2018; Lijoi, Muliere, Prünster, and Taddei, 2016). Gamma distribution also serves as the conjugate prior for many model parameters, such as a normal precision and Poisson intensity. The normalizing constant depends on gamma function whose argument is the shape parameter. Bayesian Inference in gamma models is a long standing problem that presents significant technical and computational difficulties (Damsleth, 1975; Rossell, 2009; Miller, 2018). Similar issue also occurs to the learning of other widely-used gamma models. Table 1 gives a list of distributions, where the shape/concentration/dispersion parameters are nested in gamma functions.

By exploiting normal variance-mean mixture identities related to gamma function, we derive a general data augmentation strategy. Our main result is given in Proposition 1, according to which a MCMC algorithm is built in Proposition 2 and an expectation-maximization algorithm is developed in Section 3.1.

The main difficulty is to address the reciprocal gamma function $1/\Gamma(\alpha)$. Following Barndorff-Nielsen, Kent, and Sørensen (1982), Hartman (1976) and Roynette and Yor (2005), we represent $1/\Gamma(\alpha)$ as a mean-variance mixture of normals, then we define our new class of distributions named the Pólya Inverse Gamma (P-IG)111It was previously circulated under the name of Exponential Reciprocal Gamma class. This novel distribution places gamma models in the same footing as other commonly used Bayesian models, such as sparsity (lasso, horseshoe) and logit (Pólya-Gamma). As a by-product, we show how to nest the latter within our framework. Thus we unify the inference procedure for many models.

Our data augmentation strategy with P-IG auxiliary variables may be utilized to design efficient Markov chain Monte Carlo (MCMC) algorithms in latent Dirichlet allocation (Blei, Ng, and Jordan, 2003), Beta-negative binomial models (Zhou, Hannah, Dunson, and Carin, 2012), and Gamma-Gamma (GaGa) hierarchical models (Rossell, 2009). It adds to the literature on Bayesian computation with auxiliary variables, which have proven useful in computing posterior distributions in logistic regression (Polson, Scott, and Windle, 2013), negative binomial regression (Pillow and Scott, 2012), multinomial factor models (Holmes and Held, 2006), support vector machines (Mallick, Ghosh, and Ghosh, 2005; Polson and Scott, 2011), and dependent multinomial models (Linderman, Johnson, and Adams, 2015).

To illustrate our methodology, we show examples including gamma shape inference, negative binomial regression and multinomial-Dirichlet model. The first example has a posterior density which exactly matches with the above form. Hence our algorithms can be straightforwardly applied. The posterior samples are efficiently drawn and the posterior mode is easily found by EM algorithm. The second example extends our results to incorporate Pólya-Gamma mixture representation described in (Polson, Scott, and Windle, 2013). The last example demonstrates how our methodology generalizes to high-dimensional case. Conditioned on the auxiliary variables, all elements of the multivariate concentration parameter become mutually independent, which reduces the sampling difficulty significantly. We also suggest the use of normal approximation to speed up the sampling procedure in this example.

The P-IG family of distributions is defined as an infinite convolution of Generalized inverse Gaussian (GIG) distributions and is related to the class of Pólya-Gamma (PG) distributions (Polson, Scott, and Windle, 2013) for logistic regression. Other mixture representations related to GIG are introduced in Zhang, Wang, Liu, Jordan, and Lawrence (2012) and Barndorff-Nielsen and Shephard (2012), with applications in sparse regression and stochastic volatility modelling. The P-IG(0) distribution is also a special case of $H_{a}^{\Gamma}$ distribution family studied in Roynette and Yor (2005), who provide a representation of the ratio gamma functions as a scale mixture of normals. This adds to scale mixtures results in Bayesian inference, see Andrews and Mallows (1974), Barndorff-Nielsen, Kent, and Sørensen (1982), West (1987), and Polson, Scott, and Windle (2013). Scale mixtures of normals are increasingly used in modeling complex high-dimensional distributions, and Bhattacharya, Chakraborty, and Mallick (2016) provide fast sampling strategies, adding to the practical use of scale mixture distributions in scalable stochastic simulations. Equivalently constructed scalable PG sampling schemes are provided in Windle, Polson, and Scott (2014) and Glynn, Tokdar, Howard, and Banks (2019).

1.1 Connections with Existing Work

Obtaining random draws or finding the mode from a posterior distribution involving gamma functions is computationally challenging it requires accurate gamma function evaluations. Approximation methods have been proposed to handle the computational burden. Minka (2000) describes an efficient iterative schemes for maximum likelihood estimate of Dirichlet distribution. The likelihood of Dirichlet precision is approximated by a simpler function (gamma density) by matching the first two derivatives, while the multivariate Dirichlet mean is estimated separately by fixed-point iteration. Miller (2018) applies the same idea of derivative-matching and approximate the full conditional distribution of gamma shape parameter by a gamma density function. Rossell (2009) defines a gamma shape distribution in differential expression analysis. By approximating gamma function with Stirling’s formula and evaluating the limit of the expression, the proposed distribution is roughly proportional to a gamma density. However, These ad-hoc methods are all essentially derived in univariate case. It’s not straightforward to generalize them in multivariate cases as we need to deal with correlations, and the computation of approximating parameters itself will get cumbersome. Our framework instead provides an easy and unified way to derive both MCMC and EM algorithms for multivariate models, without the need to approximate density functions.

The rest of our paper proceeds as follows: Section 2 defines the class of P-IG distributions; Section 3 illustrates our data augmentation strategy, developing a parameter expanded Gibbs sampler as well as EM algorithm; Section 4 presents examples of gamma shape inference, negative binomial regression and multinomial-Dirichlet model; and Section 5 concludes with directions for future research.

2 Pólya Inverse Gamma (P-IG) Distribution

In this section, we present the theoretical development of the P-IG distribution class, defining the P-IG distribution by the form of its integral transform. In Section 2.1, we define a baseline case of the P-IG distribution and prove that it is an infinite convolution of independent GIG distributions; in Section 2.2, the general class of exponential reciprocal distributions is constructed with an exponential tilting of the baseline case defined in Section 2.1. Other convolutions of GIG distribution are discussed in Section 2.3.

2.1 P-IG(0) Distribution

Let P-IG$(0)$ denote the Pólya Inverse Gamma distribution. The parameter, which is a tilting parameter fixed at zero in this case, will be discussed in greater detail in Section 2.2.

Definition 2.1**.**

Random variable $X_{0}$ has an Pólya Inverse Gamma distribution, P-IG$(0)$, if and only if its density function $p_{0}(x)$ satisfies the following identity

[TABLE]

where $\gamma=-\psi(1)\approx 0.57721$ is the Euler-Mascheroni constant and $\psi(s)=\frac{d}{ds}\log\Gamma(s)$ is the digamma function.

Remark 1**.**

The product representation for the reciprocal gamma function due to Weierstrass is,

[TABLE]

Remark 2**.**

Roynette and Yor (2005)** prove the existence of an infinitely divisible distribution $H_{a}^{\Gamma}$, with density function $p_{a}^{\Gamma}(x)$, such that for $a>0$

[TABLE]

which follows from the general product representation for the reciprocal gamma function,

[TABLE]

Note that the representation in Remark 1 is a special case with $a=1$. Hence

[TABLE]

Remark 3**.**

The P-IG($n,0$) is defines as follows. If $X_{n,0}\sim\text{P-IG}(n,0)$,

[TABLE]

Note that $\text{P-IG}(n,0)$ is the equivalent in distribution to the sum of $n$ independent $\text{P-IG}(1,0)$ when $n$ is a positive integer.

2.2 General P-IG(c) Distribution

We now construct the general class of P-IG distributions, P-IG$(c)$, by exponentially tilting the P-IG$(0)$ distribution. The exponential tilting strategy – similar to the one used by Polson, Scott, and Windle (2013) – allows a second parameter $c\in\mathbb{R}^{+}$ to inform a priori the precision of the P-IG random variable.

Definition 2.2**.**

The P-IG$(c)$ distribution is constructed as an exponential tilting of the P-IG$(0)$ density. Its density function is

[TABLE]

The normalizing constant, namely $Z_{c}=1/E\left(\exp\left(-c^{2}X_{0}\right)\right)$ where $X_{0}$ is an P-IG$(0)$ random variable, can be calculated using the identity in Remark 1. Similarly, the integral identity of P-IG$(c)$ is given by

[TABLE]

Our main result, presented in Theorem 1, is that a random variable $X_{c}\sim$ P-IG$(c)$ may be constructed from an infinite sum of independent generalized inverse Gaussian (GIG) random variables. The power of the result lies in the ability to identify previously unknown conditional posterior distributions.

Theorem 1**.**

The P-IG$(c)$ class of distributions can be constructed as an infinite sum of independent generalized inverse Gaussian (GIG) distributions as follows

[TABLE]

In particular, when $c=0$, the GIG distribution reduces to inverse gamma distribution. Hence,

[TABLE]

where $\Gamma^{-1}_{k}$ are i.i.d. inverse gamma random variables with shape $\frac{3}{2}$ and scale 1.

Proof.

See Appendix A. ∎

The following theorem concerns the first two moments of P-IG($c$) which will be used to construct our EM algorithm in Section 3.1 and the approximate Gibbs sampler in Appendix D.

Theorem 2**.**

If $G_{k}$ is a GIG random variable with $p=-\frac{3}{2},a=2c^{2},b=\frac{1}{2k^{2}}$, then the mean and variance of the tail infinite sum $\sum_{k=N}^{\infty}G_{k}$ are

[TABLE]

where $\psi(s)=\frac{d}{ds}\log\Gamma(s)$ is the digamma function. Setting $N=1$ gives us the first two moments of P-IG($c$).

Proof.

If $g_{k}$ is a GIG random variable with $p=-\frac{3}{2},a=2c^{2},b=\frac{1}{2k^{2}}$, then its mean and variance are given by

[TABLE]

Theorem 2 is thus a direct application of Theorem 1. ∎

Remark 4**.**

The P-IG$(c)$ distribution class belongs to the family of generalized gamma convolutions (GGC), Bondesson (1992). Its Laplace transform in Equation (4) also satisfies

[TABLE]

Its Lévy density $\nu(x)$ and Thorin density $\mu(t)$ are

[TABLE]

This result can be used to generate P-IG random variables as it shows that it falls into the class of Generalized Gamma Convolution, see Bondesson (1982) and Rosiński (2001).

2.3 GIG Mixtures

The P-IG distribution allows us to represent the unnormalized density $\frac{e^{ax}}{\Gamma(1+x)}$ as a normal variance-mean mixture. That is,

[TABLE]

where $\phi(\cdot\mid\mu,\sigma^{2})$ is the normal density with mean $\mu$ and variance $\sigma^{2}$ and $P_{0}(\cdot)$ is the distribution function of P-IG(0). $\mu(\omega)=(a+\gamma)/(2\omega),\sigma^{2}(\omega)=1/(2\omega)$ and $\tau(\omega)=\sqrt{\pi/\omega}\exp\left\{(a+\gamma)/4\omega\right\}$. In statistics and machine learning, probability distributions whose density function $p(x)$ is of the following form are used explicitly and implicitly:

[TABLE]

Here $f(x\mid\theta(\omega))$ is some well-known density function, e.g. normal, and the mixing $p(\omega)$ is the distribution of a single GIG or an infinite convolution of GIG’s. Combined with a data-augmentation scheme, the above mixture representation provides a powerful framework to solve many non-Gaussian models.

When $f=\phi$ and the mixing distribution is GIG, Polson and Scott (2013) give the variance-mean mixture representations for many common loss functions in regression and binary classification problems, which corresponds to different choices of the function $\theta(\omega)=(\mu(\omega),\sigma^{2}(\omega))$ and parameters of GIG. For example, absolute loss $L(y)=|y|$, hinge loss for support vector machine $L(y)=\max(1-y,0)$ and check loss for quantile regression $L(y)=|y|+(2q-1)y$. The representations then help reduce those non-Gaussian models to Gaussian linear models with heteroscedastic errors. Note that GIG is a very general family with many common distributions as its special cases, such as gamma, inverse gamma and inverse Gaussian distribution. 2. 2.

For the logistic loss in binary classification $L(y)=\log(1+e^{y})$, Polson and Scott (2013) show that it is also a normal variance-mean mixture. The mixing distribution is Pólya distribution, which is constructed as an infinite sum of exponentials. Note that exponential distribution is again a special case of GIG. 3. 3.

By choosing $f$ to be the exponential power density, $f(x\mid\eta(\omega),q)\propto\exp\left\{-\frac{1}{2\eta(\omega)}|x|^{q}\right\}$, and the mixing distribution $p(\omega)$ to be GIG, Zhang, Wang, Liu, Jordan, and Lawrence (2012) introduce a sparsity-inducing prior called EP-GIG and develop EM algorithms for sparse learning. The density function of EP-GIG is given explicitly and special cases (generalized $t$ distribution and exponential power-gamma distribution) are discussed when the mixing GIG reduces to inverse gamma and gamma respectively. 4. 4.

The Pólya-Gamma distribution class is proposed by Polson, Scott, and Windle (2013) to solve the inference problem in models with binomial likelihoods, including logistic regression and negative binomial regression. Pólya-Gamma distribution PG$(b,c)$ can be written as an infinite convolution of gamma distributions whose shapes are all equal to $b$ and scales depend on $c$, or equivalently $\text{GIG}(b,2,0)$.

[TABLE]

Choose $f(z\mid\omega)=2^{-b}\exp\left\{-\omega z^{2}/2+\kappa z\right\}$ with $\kappa=a-b/2$ and $z=x^{\prime}\beta$, the likelihood of logistic regression is a mixture

[TABLE]

Here $p_{\text{PG}}(\omega\mid b,0)$ is the density of PG$(b,0)$ and $f(z\mid\omega)$ is proportional to the normal density where the mean and variance are functions of $\omega$. 5. 5.

Barndorff-Nielsen and Shephard (2012) use a normal variance-mean mixture as a general approach of building densities on the real line. Here $f(x\mid\theta(\omega))=\phi(x\mid\mu+\beta\omega,\omega)$. When $p(\omega)$ is GIG density, the resulted mixture is generalized hyperbolic distribution, which includes many special cases such as normal inverse Gaussian, normal gamma, Laplace, skewed Student’s $t$ distribution. Furthermore, normal distribution can also be written as a limiting case of generalized hyperbolic distribution.

3 MCMC and Data Augmentation

This section illustrates the data augmentation strategy and sampling scheme for Gamma inference using the P-IG class distributions. First, notice that many Bayesian gamma models involve a posterior density of the form

[TABLE]

where $x\in\mathbb{R}^{+}$ is on the positive real line, and $C$ is the normalizing constant. The arguments in gamma functions, $\left\{g_{\ell}(\cdot)\right\}_{\ell=1}^{L},\left\{h_{m}(\cdot)\right\}_{m=1}^{M}$ and $\left\{j_{n}(\cdot)\right\}_{n=1}^{N}$ are nonnegative increasing linear functions of $x$ on $(0,\infty)$. We assume $M\geq 1$, otherwise the form might not be integrable. Parameters $p,a,b,\beta_{1},...,\beta_{N}$ are scalars. $\beta_{n}$’s are positive. $p$ and $a$ are nonnegative. In order to perform full posterior inference on the variable $x$, a sampling procedure for $p(x)$ is needed, while a Maximum A Posteriori (MAP) estimate requires finding the maximizer of it.

The idea of the data augmentation strategy is to introduce a group of auxiliary random variables $\bm{\omega}=(\omega_{1},...,\omega_{n})^{T}$ such that for $\bm{\omega}\in\Omega$,

[TABLE]

and the joint density $p(x,\bm{\omega})$ after augmentation is easier to deal with, as it doesn’t consist of gamma functions any longer.

Returning to Equation (7), the posterior density involves gamma functions $\Gamma(\cdot)$, reciprocal gamma functions $1/\Gamma(\cdot)$ and gamma ratios $\Gamma(\cdot)/\Gamma(\cdot+\beta)$. We will then express each of them using the corresponding integral representation in Table 2. This is equivalent to introducing an auxiliary random variable for each of them. The total number of the auxiliary random variables is then $L+M+N$.

Proposition 1**.**

The posterior density admits an integral representation as follows

[TABLE]

where

[TABLE]

$p_{0}(\cdot)$ is the probability density of P-IG(0).

Remark 5**.**

When $N\geq 1$ in the form (7), the third term $\left(\prod_{n=1}^{N}\frac{\Gamma(j_{n}(x))}{\Gamma(j_{n}(x)+\beta_{n})}\right)$ can be absorbed into the first two terms. However, we still recommend using the Beta representation if possible. Otherwise the total number of auxiliary variables needed is increased by $N$.

Remark 6**.**

We may generalize the form (7) in a few ways:

For multivariate $x$ of dimension $d$, if $g,h,j$ are all linear functions mapping from $\mathbb{R}^{+,d}$ to $\mathbb{R}^{+}$, then the data augmentation is still valid. 2. 2.

The $x^{p}$ term can be replaced with a polynomial function of $x$, as long as it’s always positive for $x>0$. 3. 3.

If the posterior density $p(x)$ has extra factors which are not included in the form (7), the strategy still works as long as we can find integral representations for those extra factors.

Although the posterior joint distribution in Equation (8) looks forbidding at the first glance, in many applications the linear functions $g,h,j$ are simple enough, e.g. $h_{m}(x)=x$ for all $m$, and $N=0$, which simplifies the expression a lot. More importantly, the conditional posteriors can be derived easily. To derive that of $\tau_{\ell}$, for example,

[TABLE]

notice that it’s proportional to $\tau_{\ell}^{g_{\ell}(x)-1}e^{-\tau_{\ell}}$. Therefore the conditional posterior distribution of $\tau_{\ell}$ is $\Gamma(g_{\ell}(x),1)$. Similarly for $\omega_{m}$’s and $\eta_{n}$’s. For the conditional posterior of $x$, since $g,h,j$ are all linear, it’s proportional to $Q_{1}(x)e^{-Q_{2}(x)}$ where $Q_{1}(x)$ is a polynomial of degree $p+M-1$ and $Q_{2}(x)$ is a quadratic function.

Before we summarize the conditional posteriors with respect to the augmented vector $(x,\bm{\tau,\omega,\eta})$, the following definition of power truncated normal (PTN) distribution is useful.

Definition 3.1**.**

The power truncated normal distribution, PTN$(p,a,b)$, has density function

[TABLE]

where $p,a>0$ and $b\neq 0$.

Finally, the following proposition is helpful in developing the Gibbs sampler.

Proposition 2** (Gibbs Sampler).**

If the joint probability density of the augmented vector $(x,\bm{\tau,\omega,\eta})$ is proportional to $G(x,\bm{\tau})\cdot H(x,\bm{\omega})\cdot J(x,\bm{\eta})\cdot Q(x)$ given in Proposition 1, then the conditional distributions are

[TABLE]

Here the conditional distribution of $x$ is a finite discrete mixture of PTN distributions. $\{\pi_{m}\}_{m=0}^{M}$ are proportional to the coefficients in the polynomial $Q(x)$.

[TABLE]

Furthermore, given $x$, all auxiliary random variables are conditionally independent. The sampling methods for P-IG and PTN are given in Appendix B and C.

In many statistical applications, the number of auxiliary random variables grows linearly with the sample size and problem dimension, but the forms of $g,h,j$ are relatively simple. Observing that $\tilde{a}$ and $\tilde{b}$ are both the sum of numerous terms, we may use normal variables to approximate them, by matching the moments, so that the sampling procedure for $(\bm{\tau,\omega,\eta})$ can be skipped in Gibbs sampling. Appendix D gives the approximate Gibbs sampler when $g_{\ell}(x)=g(x)$ for all $\ell$ (similarly for $h$ and $j$), and $L,M,N$ are all large enough.

3.1 Expectation-Maximization Algorithm

By exploiting Proposition 1, MAP of the model can be found using an expectation-maximization algorithm. The complete-data log posterior is

[TABLE]

for some constants $c_{0},c_{1}$ (with respect to $x$). In the $t$-th expectation step, we compute the expected value of $\log p(x,\bm{\tau,\omega,\eta})$ under the current conditional posterior $p(\bm{\tau,\omega,\eta}\mid x_{(t)})$, denoted as $C(x\mid x_{(t)})$. Then in the maximization step, $C(x\mid x^{(t)})$ is maximized as a function of $x$. We now derive the expectation and maximization steps.

• •

**The Expectation Step

**From equation (10), notice that $\log p(x,\bm{\tau,\omega,\eta})$ is linear in terms of $\log\tau_{\ell},\omega_{m}$ and $\log\eta_{n}$. Therefore, we replace them with their conditional expectations in the expectation step. Applying Proposition 2 yields

[TABLE]

The derivation above uses the properties of gamma and beta distribution, as well as Theorem 2 which calculates the expectation of P-IG($c$) distribution. Finally, one can represent the function $C(x\mid x_{(t)})$ (up to a constant) as

[TABLE]

Given the linearity of $g,h,j$ functions, it can be further simplified as

[TABLE]

for some constant $\kappa_{1}>0$ and $\kappa_{2}\in\mathbb{R}$, depending on $x_{(t)}$.

• •

**The Maximization Step

**Since $C(x\mid x_{(t)})\rightarrow-\infty$ as $x\rightarrow\infty$, we conclude that $C(x\mid x_{(t)})$ as a function of $x$ has a maximizer on $(0,\infty)$, which can be found numerically. Furthermore, when $h_{m}(0)=0$ for all $m$, the unique maximizer has a closed form

[TABLE]

4 Examples

4.1 Inference for Gamma Shape

The gamma distribution, parameterized by shape $\alpha$ and rate $\beta$, is a component of many probability models. For instance, a gamma prior distribution for the precision parameter in Gaussian linear models is quite common. In fact, Normal-gamma distributions are workhorse models for shrinkage estimation in regression problems (Griffin and Brown, 2010). Gamma distribution is also widely used in modelling of extreme values, where it serves as the prior for the shape parameter of Pareto distribution (Arnold and Press (1989)) and helps to construct a quasi-conjugate prior for generalized Pareto distribution (Diebolt, El-Aroui, Garrido, and Girard (2005)). While a gamma prior distribution for a parameter is common, it is less common to model hyperparameters of the gamma distribution itself as random variables – particularly the shape parameter, $\alpha$. Although posterior inference of the rate parameter is straightforward – since the gamma distribution itself is a conjugate prior for the rate parameter – posterior inference of the gamma shape parameter is a long-standing problem (Damsleth, 1975; Damien, Laud, and Smith, 1995; Rossell, 2009; Miller, 2018) and efficient posterior computation remains an open problem.

Damsleth (1975) discussed two conjugate priors for $\alpha$, the gamma shape parameter. They are called GamCon distributions of type I and type II. Type I assumes the rate $\beta$ is known, while the type II doesn’t. In this section, we focus on the latter one and replicate Damsleth’s example, showing how to utilize the data augmentation scheme and build algorithms for posterior inference.

Let’s consider the following hierarchical model

[TABLE]

Here $\xi_{2}$ is GamCon distribution of type II. $\eta>\mu>0,\delta>0$. The probability density of $\xi_{2}(\mu,\delta)$ is

[TABLE]

Suppose observations $\bm{x}=(x_{1},x_{2},...,x_{n})^{\prime}$ are independently and identically distributed gamma random variables drawn from the above model. The likelihood is

[TABLE]

where $x_{g}$ is the geometrical mean, $x_{g}=\left(\prod_{i=1}^{n}x_{i}\right)^{1/n}$.

The marginal posterior distribution of $\alpha$ is then calculated as

[TABLE]

By construction, the marginal posterior of $\alpha$ given $\bm{x}$ also follows $\xi_{2}$, with updated parameters $(\eta^{\prime},\mu^{\prime},\delta^{\prime})$. That is, $\alpha\mid\bm{x}\sim\xi_{2}(\eta^{\prime}/\mu^{\prime},\delta^{\prime})$ and

[TABLE]

where $x_{a}$ is the arithmetical mean, $x_{a}=\frac{1}{n}\sum_{i=1}^{n}x_{i}$. One can observe that $\eta^{\prime}$ is a weighted arithmetical mean of $\eta$ and $x_{a}$ with weights $\delta$ and $n$ respectively and $\mu^{\prime}$ is a weighted geometrical mean of $\mu$ and $x_{g}$, also with weights $\delta$ and $n$. Here $x_{a}$ and $x_{g}$ are two sufficient statistics. $\delta$ is viewed as the prior sample size while $\eta$ and $\mu$ are the prior means. The problem of gamma shape inference is thus translated to that of $\xi_{2}$ distribution.

4.1.1 MCMC and EM for $\xi_{2}$ when $\delta\in\mathbb{N}^{+}$

We first develop the Gibbs sampler for a general GamCon distribution of type II, using data augmentation strategy. Suppose $x\sim\xi_{2}(x;\mu,\delta)$ and $\delta$ is a positive integer. Then the probability density can be rewritten in the form (7)

[TABLE]

Hence $L=1$ with $g_{1}(x)=\delta x+1$, $M=\delta$ with $h_{m}(x)=x$ and $N=0$. $p=1,a=0$ and $b=-\delta\log(\delta\mu)$.

After introducing the auxiliary variables $\tau$ and $\omega_{1},\omega_{2},...,\omega_{\delta}$, Proposition 2 then immediately gives the conditional posteriors:

[TABLE]

4.1.2 Damsleth Examples with Non-Informative Prior

Here we replicate Damsleth’s example of no prior information with sample size $n=5,10,30$. Putting $\delta=0$, the posterior parameters are

[TABLE]

The resulted posterior of gamma shape is thus $\alpha\mid\bm{x}\sim\xi_{2}(x_{a}/x_{g},n)$. Instead of actually generating Gamma random variables, we directly use the sufficient statistics $(x_{a},x_{g})$, given in Table 2 of Damsleth (1975). The true value of $\alpha$ is 5. The histogram of 5000 posterior samples for each case are shown in Figure 1. The colored solid lines denote the true posterior density and the black dashed line denotes the true $\alpha$. We see that the posterior samples are indeed sampled from the target distributions. Figure 2 shows the sampling trace plots. In Table 3, the sample moments are compared with their theoretical counterparts calculated by numerical integration. When $n=30$, the deviations from the theoretical values are $-1.0\%,-4.6\%,-8.2\%$ and $0.1\%$ for mean, variance, skewness and kurtosis respectively.

To find the posterior mode of $\alpha$, we exploit the EM algorithm developed in Section 3.1. Specially, in the expectation step, the conditional expected value of complete-data log posterior (at $t$-th iteration) given by Equation (11) now has the parameters

[TABLE]

Then in the maximization step, $\alpha$ is updated by Equation (12), with $p-M+1=\delta^{\prime}$. In Figure 3, we start with 30 different initial values in EM algorithm and show the optimizing paths for the case $n=30$. In this particular case, the algorithm convP-IGes after around 50 steps. The 30 numerical solutions given by EM has mean 3.054, which matches with the maximizer found by Mathematica. And the standard deviation is as small as 0.0003.

4.2 Negative Binomial Regression

Next, we proceed to how our data augmentation strategy can be used to fit negative binomial regression models. The count data $\{y_{i}\}_{i=1}^{n}$ are assumed to follow the negative binomial distribution

[TABLE]

where $\{{\bm{x}_{i}}\}_{i=1}^{n}$ are observed covariates and $\bm{\beta}$ is the regression coefficients. $r$ is interpreted as the number of failures until the experiment is stopped, while the success probability $p_{i}$ is related to $\bm{x_{i}^{\prime}\beta}$ via the logistic transformation. This model specification is equivalent to the following Gamma-Poisson mixture,

[TABLE]

The likelihood is

[TABLE]

Pillow and Scott (2012) adopt exactly the same model with known $r$ and use a data augmentation strategy for the inference on $\beta$. Our example here extends their method and allows for inference on both $r$ and $\beta$. The parameter $r$ controls the dispersion of observations, as the expectation of $y_{i}$ is $re^{\bm{x_{i}^{\prime}\beta}}$ and variance is $re^{\bm{x_{i}^{\prime}\beta}}\left(1+e^{\bm{x_{i}^{\prime}\beta}}\right)$.

Let the prior for $\bm{\beta}$ and $r$ be $N(\bm{0},\Sigma)$ and $p(r)$. We then calculate the joint posterior of $({\bm{\beta}},r)$ as

[TABLE]

Notice that if we assign gamma prior for $r$, then its conditional posterior density is close to the form (7), except for the extra factor $\prod_{i=1}^{n}\left(e^{\bm{x_{i}^{\prime}\beta}}\right)^{y_{i}}/\left(e^{\bm{x_{i}^{\prime}\beta}}+1\right)^{r+y_{i}}$. Similarly for $\beta$, the posterior is close to normal density, except for the same factor. This is however not an issue as we can write this factor as a scale mixture of normals, where the mixing distribution is Pólya-Gamma from Polson, Scott, and Windle (2013). Setting $a=y_{i}$ and $b=r+y_{i}$, the key mixture representation in Equation (6) related to it now becomes

[TABLE]

The integrand $f(r,\bm{\beta}\mid\xi)$ is an exponential density when viewed as a function of $r$, and a normal density when viewed as a function of $\bm{\beta}$.

Finally, we introduce the auxiliary variables $\bm{\tau,\omega,\xi}$ which follow gamma, P-IG and PG distribution respectively. Let $\Xi:=\text{diag}(\xi_{1},...,\xi_{n})$ and $\bm{z}:=\left(\frac{y_{1}-r}{2\xi_{1}},...,\frac{y_{n}-r}{2\xi_{n}}\right)^{\prime}$ and $p(r)\sim\Gamma(a_{0},b_{0})$. The conditional posteriors are derived as follows:

[TABLE]

where

[TABLE]

Figure 4 shows an illustrating simulation example where we set true $r=5$ and generate $n=100$ count observations. We draw the coefficient vector $\bm{\beta}\in\mathbb{R}^{5}$ and covariates $x_{ij}\stackrel{{\scriptstyle i.i.d.}}{{\sim}}N(0,0.5^{2})$ for $1\leq i\leq n$ and $1\leq j\leq 5$. The prior for $r$ is $p(r)\sim 1/r$ which is the limit case of gamma prior with $a_{0}=0,b_{0}=0$. For $\beta$, we choose $\Sigma=10^{6}I_{5}$ so that the prior information is relatively weak. The left panel shows the boxplots of posterior $\beta$ samples, with those red dots denote the true $\beta$’s. As the figure shows, all 5 $\beta$’s fall in the 95% credible interval. The middle panel is the histogram of the posterior samples for $r$. Since the actual marginal posterior of $r$ is hard to compute given the complicated form of the joint $p({\bm{\beta}},r\mid{\bm{X},y})$, in the right panel we use the same dataset, plug in the true $\beta$’s, and run the Gibbs sampler again, so that we are able to compare the resulted histogram with the true density (red line). They are quite close to each other, indicating that the sampling procedure works well for this example.

4.3 Multinomial-Dirichlet Model

In this section, we develop the Markov chain Monte Carlo (MCMC) algorithm for fully posterior inference of the concentration parameter vector in the Dirichlet distribution. Such inference problems commonly arise in applied analyses of categorical data. Section 4.3.1 presents the general hierarchical multinomial-Dirichlet model class for which the P-IG data augmentation scheme may be utilized. Section 4.3.2 develops a Gibbs sampler for inferring the concentration parameter in the Dirichlet distribution and conducts a simulation study comparing the performance of our data augmentation strategy with Metropolis-Hasting algorithm.

4.3.1 A Hierarchical Multinomial-Dirichlet Model

The multinomial-Dirichlet framework presented herein is closely related to the latent Dirichlet allocation model of Blei, Ng, and Jordan (2003) for topic modeling of text data, and we use text analysis as a motivating context. Suppose that for document $s\in\{1,...,S\}$, each of $N_{s}$ words in the document is independently allocated to $K$ topics conditional on probability vector $\bm{p}_{s}=\left(p_{s1},p_{s2},\ldots,p_{sK}\right)$. For each document $s$, the number of words allocated to each topic is denoted by $\bm{n}_{s}=\left(n_{s1},\cdots,n_{sK}\right)$, which follows a multinomial distribution. The sampling model for the count vector $\bm{n}_{s}$ is then a multinomial distribution given probability vector $\bm{p}_{s}$,

[TABLE]

The probability vector $\bm{p}_{s}$ is the proportional allocation of each document to the $K$ topics. In a Bayesian analysis, the probability vector for each document $\bm{p}_{s}$ is typically assigned a Dirichlet distribution with concentration parameter vector $\bm{\alpha}=\left(\alpha_{1},\ldots,\alpha_{K}\right)$,

[TABLE]

Rather than fixing $\bm{\alpha}=\left(\frac{1}{K},\ldots,\frac{1}{K}\right)$, as is common, we complete the model with a prior distribution $p(\bm{\alpha})$. This hierarchical prior distribution for $\bm{\alpha}$ facilitates more efficient information sharing across documents (observational units), and it yields practical advantages for out-of-sample prediction, which we discuss below. The model framework and P-IG augmentation admit independent uniform, truncated normal, and exponential prior distributions for the elements $\alpha_{k}$. Section 4.3.2 presents analyses based on independent gamma priors $p(\bm{\alpha})=\prod_{k=1}^{K}p(\alpha_{k})$.

In application, model inferences are often summarized by the posterior predictive distribution for the topic proportion vector $\bm{p}^{*}$ in a new document. Computing the posterior predictive distribution $p\left(\bm{p}^{*}\mid\bm{n}_{1},...,\bm{n}_{S}\right)=\int_{\bm{\alpha}}p\left(\bm{p}^{*}\mid\bm{\alpha}\right)p(\bm{\alpha}\mid\bm{n}_{1},...,\bm{n}_{S})d\bm{\alpha}$ requires posterior computation of $p\left(\bm{\alpha}\mid\bm{n}_{1},\ldots,\bm{n}_{S}\right)\propto p\left(\bm{\alpha}\right)\prod_{s=1}^{S}p\left(\bm{n}_{s}\mid\bm{\alpha}\right)$; however, when the probability vectors $\bm{p}_{s}$ are integrated out of the multinomial likelihood, the marginal likelihood $p\left(\bm{n}_{s}\mid\bm{\alpha}\right)$ includes elements of $\bm{\alpha}$ inside the gamma function,

[TABLE]

Because $\bm{\alpha}$ is nested inside the gamma function, computing $p(\bm{\alpha}\mid\bm{n}_{1},\ldots,\bm{n}_{S})$ is a challenge. Previous inference strategies relied on approximations, but in Section 4.3.2 we introduce a new data augmentation scheme for computing the full posterior $p\left(\bm{\alpha}\mid\bm{n}_{1},\ldots,\bm{n}_{S}\right)$.

4.3.2 Data Augmentation and Simulation Study

Assume independent gamma prior distributions for each element of vector $\bm{\alpha}$ so that $p\left(\bm{\alpha}\right)\propto\prod_{k=1}^{K}\alpha_{k}^{a_{0}-1}e^{-b_{0}\alpha_{k}}$, with hyperparamter $a_{0}$ and $b_{0}$. Note that gamma priors on each $\alpha_{k}$ give closed-form full conditional distributions in a Gibbs sampler, which is shown below. When $\alpha_{k}\sim\Gamma(a_{0},b_{0})$, where $a_{0}$ denotes the shape parameter and $b_{0}$ the rate, the expectation $E[\alpha_{k}]=a_{0}/b_{0}$. We can set $E[\alpha_{k}]=1/K$, a standard choice for the Dirichlet concentration parameter, by choosing $a_{0}=b_{0}/K$. The prior variance then depends on both the dimension of the Dirichlet distribution, $K$, and the rate parameter, $b_{0}$.

We now reorganize the multivariate posterior density of $\bm{\alpha}$ as

[TABLE]

where

[TABLE]

Note that for each $\alpha_{k}$, Equation (15) is exactly of the form (7). And the extra factor, $\prod_{s=1}^{S}\frac{\Gamma(\sum_{k}\alpha_{k})}{\Gamma(\sum_{k}\alpha_{k}+N_{s})}$, can be replaced with a beta integral representation. A multivariate version of Proposition 2 produces the conditional posteriors as follows

[TABLE]

where $a_{k}=\sum_{s=1}^{S}\omega_{sk},b_{k}=S\gamma-b_{0}+\sum_{s=1}^{S}\log(\tau_{sk}\cdot\eta_{s})$. Conditioned on the introduced auxiliary variables $(\bm{\tau,\omega,\eta})$, all elements of the vector $\bm{\alpha}$ are now mutually independent, which significantly reduces the difficulty of sampling procedure as we can now sample separately from $K$ univariate distributions.

The total number of auxiliary variables to be sampled at each iterate is $S(1+2K)$, which may greatly slow down the MCMC algorithm for large values of $S$ and $K$. However, we observe that these variables affect the distribution of $\alpha_{k}$ only through the parameters $(a_{k},b_{k})$. The summation form of $(a_{k},b_{k})$ suggests that we may apply central limit theorem and approximate them by normal variables. Therefore, as mentioned in Section 3, the Gibbs sampling algorithm can be approximately simplified to

Initialize $a_{k}^{(0)},b_{k}^{(0)}$ for $1\leq k\leq K$. 2. 2.

At step $t$, sample $a_{k}^{(t)}$ and $b_{k}^{(t)}$ from $N(\mu_{a,k},\sigma^{2}_{a,k})$ and $N(\mu_{b,k},\sigma^{2}_{b,k})$ respectively, where

[TABLE]

and $\alpha_{0}=\sum_{k=1}^{K}\alpha_{k}$. 3. 3.

Sample $\alpha_{k}^{(t)}$ from $\text{PTN}(S+a_{0},a_{k}^{(t)},b_{k}^{(t)})$ for $1\leq k\leq K$. Increase $t$ by 1 and return to (2).

Fixing the dataset dimensions $S$ and $K$, we consider two settings of the true $\bm{\alpha}$ for our simulation experiment: (A) heterogeneous $\alpha_{k}=k/K$ and (B) homogeneous $\alpha_{k}=1/K$ for $k=1,2,...,K$. Probability vector $\bm{p}_{s}$ are drawn independently from $\text{Dirichlet}(\bm{\alpha})$ and vector of counts $\bm{n}_{s}$ are drawn from $\text{Mutinomial }(N_{s},\bm{p}_{s})$ with $N_{s}=500$ for $s=1,2,...,S$. We set the hyperparamters $a_{0}=b_{0}/K$ in independent gamma priors of $\alpha_{k}$’s, so that the prior expectations are all equal to $1/K$. The setting of our simulation experiment is summarized below:

• •

$S\in\{100,1000\},K\in\{10,50\}$ and $b_{0}\in\{0.1,1,5\}$.

• •

4 MCMC algorithms in comparison:

• –

DA: Gibbs sampler which iteratively samples $\bm{\tau,\omega,\eta}$ and $\bm{\alpha}$.

• –

DA-N: replace the parameters in the conditional posterior of $\bm{\alpha}$ with approximating normal variables and skip the sampling of $\bm{\tau,\omega,\eta}$ in DA.

• –

DA-E: replace the parameters in the conditional posterior of $\bm{\alpha}$ with corresponding expectations and skip the sampling of $\bm{\tau,\omega,\eta}$ in DA.

• –

MH: random-walk Metropolis-Hasting sampler.

• •

Metrics:

• –

Root Mean Square Error:

[TABLE]

$T=500$ is the posterior sample size.

• –

Effective sample size ratio (ESSR):

[TABLE]

which measures the serial correlation between posterior samples. $\lambda^{2}$ is the sample variance and $\sigma^{2}$ is the estimate of the spectral density at frequency zero.

• –

Algorithm running time in seconds

In simulation experiments, we initialize each $\alpha_{k}$ with a random draw from $\text{Lognormal}(0,1/K)$ and multiply by $1/K$, so that the prior median is $1/K$. For each combination of $(b_{0},S,K)$, we run the simulation for 50 times (the first 200 samples are dropped each time). Table 4 and 5 show the results, averaged over 50 runs, for the homogeneous and heterogeneous setting respectively. Despite being slow, our data augmentation strategy with P-IG auxiliary variables produces more accurate estimates of $\bm{\alpha}$ than Metropolis-Hasting does, in heterogeneous setting. While in homogeneous setting, the two methods have similar RMSE. Metropolis-Hasting gets significantly worse when $K$ grows to 50. This is not surprising at all. As we notice that, the effective sample size ratios for Metropolis-Hasting are as low as 0.01, indicating strong serial correlations in posterior samples drawn by Metropolis-Hasting, of which the acceptance rate is around 0.02. Therefore, the RMSE for Metropolis-Hasting is entirely up to the distance between the initial $\bm{\alpha}_{(0)}$ and the true $\bm{\alpha}$. Our Gibbs sampler instead enjoys the advantage of being less sensitive to initialization and not requiring further tuning. The different choices of $b_{0}$ don’t seem to affect the posteriors too much. Once replacing auxiliary variables with normal approximations or expectations, the Gibbs sampler gets much faster. Meanwhile, root mean squared errors and effective sample sizes get slightly improved. With DA-N and DA-E, we can expect the resulted posterior samples to have slightly less variation as well as less correlation.

5 Discussion

The class of Pólya Inverse Gamma (P-IG) distributions are developed as mixing distributions for models with Gamma functions. This adds to the literature on normal variance-mean mixtures by showing that they extend to a wide class of applications. Our ensuing data augmentation strategy facilitates full posterior inference for parameters in models which were hitherto hard to provide inference. The algorithms are scalable and are a fast efficient simulation method for drawing from posterior distributions with applications to many area, such as non-parametric Bayes, latent Dirichlet allocation, Gamma-Gamma hierarchical models, extreme value models, and many other Bayesian mixture models.

The focus of our paper is on theoretical and algorithmic development of P-IG auxiliary variables. Our work contributes to the literature on scale mixtures of normals (see, e.g., (Andrews and Mallows, 1974; West, 1987; Polson, Scott, and Windle, 2013)). We believe that the computational strategies developed here will provide the foundation for new and richly structured hierarchical gamma models. Applied Bayesian analyses of categorical data will benefit from increased model flexibility and information borrowing strategies.

There are a number of avenues for future research. In particular, regularized scale allocation models can be implemented using data augmentation methods of Polson and Scott (2013) with P-IG distribution. Barndorff-Nielsen, Blaesild, and Seshadri (1992) provide multivariate GIG distribution theory and relationships with Poisson processes.

Appendix A Proof of Theorem 1

The generalized inverse Gaussian distribution, $\text{GIG}(p,a,b)$, has probability density function

[TABLE]

It suffices to show that if $G_{k}\sim GIG\left(-\frac{3}{2},2c^{2},\frac{1}{2k^{2}}\right)$, then the following integral identity holds,

[TABLE]

The density of $G_{k}$ given by

[TABLE]

with normalizing constant,

[TABLE]

It follows by the algebraic calculation,

[TABLE]

as required.

Appendix B Simulating P-IG Random Variables

We consider below 3 different ways to generate independent random variables from P-IG distribution.

1. (a)

Since Theorem 1, we can approximate an P-IG$(c)$ with the finite sum

[TABLE]

where the gamma random variable $\Gamma(a_{N},b_{N})$ are to approximate the tail part by matching the first two moments given in Theorem 2. The shape and rate parameter are

[TABLE] 2. (b)

Since the generation of i.i.d. inverse gamma variables is faster than that of GIG variables, for small values of $c$, we may consider rejection sampling with P-IG(0) as the proposal density.

Generate a sample $W$ from P-IG(0) using the method in (a) and $U$ from Unif[0,1]. 2. 2.

If $U<e^{-c^{2}W}$, accept $W$ as a sample drawn from P-IG$(c)$. Otherwise, reject $W$ and return to the sampling step. 3. (c)

McLeish (2014) shows that one can simulate random variables using the saddlepoint approximation, given the cumulant generating function $k(t)=\log E(e^{tW})$ is known. The saddlepoint approximation is

[TABLE]

where $t$ solves $k^{\prime}(t)=x$. $\Phi(\cdot)$ and $\phi(\cdot)$ are the cdf and density of standard normal distribution. We can first generate a random variable $T$ using inverse transform method, such that its cdf $F(t)=\Phi(w(t))+\phi(w(t))\left(\frac{1}{w(t)}-\frac{1}{u(t)}\right)$, then $W=k^{\prime}(T)$ has cdf given by the saddlepoint approximation above.

Genearte a sample $U$ from Unif[0,1]. 2. 2.

Solve $U=F(t)$ using Newton-Raphson iteration

[TABLE] 3. 3.

$W=k^{\prime}(t)$.

Appendix C Simulating PTN Random Variables

The probability density of a random variable $\text{PTN}(p,a,b)$ is given as

[TABLE]

Note that we can write it as a multiplication

[TABLE]

where $g(x|p,b)$ is the density of a gamma random variable with shape $p$ and rate $\tau|b|-b>0$.

[TABLE]

We sample from $p(x\mid p,a,b)$ by rejection method:

Generate $X\sim\Gamma(p,\tau|b|-b)$ and $U\sim Unif[0,1]$ 2. 2.

Return $X$ until $U\leq e^{-a(X-\tau|b|/2a)^{2}}$.

If $b>0$, we set $\tau=\sqrt{1/4+2ap/b^{2}}+1/2$; if $b<0$, we set $\tau=\sqrt{1/4+2ap/b^{2}}-1/2$. Then $e^{-a(E(X)-\tau|b|/2a)^{2}}=1$.

Appendix D Approximating Gibbs Sampler

If $g_{\ell}(x)=g(x)$ for all $\ell$ (similarly for $h$ and $j$), and $L,M,N$ are all very large, then an approximate Gibbs sampler based on the conditional posterior in Proposition 2 is

Initialize $\tilde{a}^{(0)},\tilde{b}^{(0)}$. 2. 2.

At step $t$, sample $(\tilde{a}^{(t)},\tilde{b}^{(t)})$ from $N(\mu_{a}(t),\sigma^{2}_{a}(t))$ and $N(\mu_{b}(t),\sigma^{2}_{b}(t))$ respectively, where

[TABLE] 3. 3.

Sample $x^{(t)}$ from $\sum_{m=0}^{M}\pi_{k}\cdot\text{PTN}(p+m,\tilde{a}^{(t)},\tilde{b}^{(t)})$.