Bayesian mixture modeling using a mixture of finite mixtures with normalized inverse Gaussian weights

Let each observations $Y_{1},\dots,Y_{n}$ be univariate or multivariate in an Euclidean space. A MFM model is a statistical model defined by the following hierarchical representation.

where $f(y|\theta_{i})$ is a parametric model with parameter $\theta_{i}$, $\delta_{x}(\cdot)$ is a point mass at $x$ and $\xi$ is a element of a vector $\xi=(\xi_{1},\dots,\xi_{M})^{\top}$. $p_{0}(\xi)$ is a prior density function of the parameter in mixture components. The parameter space is denoted by $\Theta\subset\mathbb{R}^{d}$. Each $c_{i}$ is a latent allocation indicates to which component each observation is allocated. Given the number of components $M$, $P_{\pi}(\pi|M)$ is a probability distribution on the $M$-dimensional unit simplex, in that, this is the prior distribution for the mixing weights. We focus the case where the mixing weights can be further hierarchically expressed as follows:

where $S_{1},\dots,S_{M}$ is unnormalized weights and $\pi$ is (normalized) mixing weight vector. This representation is the most basic way to generate a probability distribution on the $d$-dimensional unit simplex $\mathbb{S}_{d}:=\{w=(w_{1},\dots,w_{d});w_{i}\geq 0,\sum_{i=1}^{d}w_{i}=1\}$. One of the most famous examples of a distribution on $\mathbb{S}_{d}$ is the Dirichlet distribution. If $S_{m}|M\sim\text{Gamma}(\gamma_{i},1)$ in (2.2), $\pi\sim\text{Dirichlet}(\gamma_{1},\dots,\gamma_{d})$, where $\text{Gamma}(a,b)$ is the gamma distribution with shape parameter $a>0$ and scale parameter $\beta>0$, and $\text{Dirichlet}(\gamma_{1},\dots,\gamma_{d})$ is the Dirichlet distribution with parameter $\gamma=(\gamma_{1},\dots,\gamma_{d})>0$. If $\gamma_{1}=\cdots=\gamma_{d}=:\gamma$, the distribution obtained by normalization is $\text{Dirichlet}(\gamma,\dots,\gamma)$. The symmetric Dirichlet distribution is often used in the framework of MFM, because of conjugacy with categorical distributions and simplicity of computation. The symmetric structure of $\text{Dirichlet}(\gamma,\dots,\gamma)$ is also essential for marginalizing out $M$ and deriving the exchangeable partition distribution. However, it is known that the estimation result is sensitive to the choice of the shape parameter $\gamma$. For example, Miller and Harrison (2018) recommended to use $\gamma=1$ as a default choice and the value works well in many cases. The use of a small $\gamma$ was recommended in terms of sparse or generalized MFMs (Rousseau and Mengersen, 2011; Malsiner-Walli et al., 2016; Frühwirth-Schnatter et al., 2021). A small value of $\gamma$ indicates that many empty components can be created. Although most previous studies allow for the appearance of the empty components, there are several problems. First, $M$ tends to be overestimated, making the interpretation of the model difficult. Second, the existing of the empty components decreases the predictive performance of the model, as we will see in a later section through density estimations. Finally, they lead to an increase in computation time. When we focus on the estimation of the number of components $M$, it is desirable to have few empty components. In other words, the discrepancy between $M$ and $k$ should be small. To achieve this, we use the normalized inverse Gaussian distribution as mixing weights, instead of the Dirichlet distribution.

We give an equivalent representation of the MFM. We can construct a discrete measure $P(\cdot)=\sum_{m=1}^{M}\pi_{m}\delta_{\xi_{m}}(\cdot)$ in the parameter space $\Theta$, almost surely, where $M$, $\xi$ and $\pi$ are realizations from the distributions $q_{M}$, $p_{0}$ and $P_{\pi}$, respectively. Let $\theta_{1},\dots,\theta_{n}$ be random variables according to $P$. Then, the model (2.1) is equivalent to the following hierarchical representation using a random measure $P$.

where $\mathcal{P}$ is a probability distribution of $P$ with parameters $q_{M}$, $h$ and $p_{0}$. The representation (2.3) is the MFM described as a nonparametric Bayesian framework with infinite mixtures. If we replace $\mathcal{P}$ with the Dirichlet process, (2.3) represents the well-known Dirichlet process mixture models (e.g., Escobar and West, 1995). If we replace $\mathcal{P}$ with the normalized inverse Gaussian process, (2.3) represents the normalized inverse Gaussian process mixture models (Lijoi et al., 2005). For MFM, Argiento and De Iorio (2022) proposed a normalized independent finite point process (Norm-IFPP), which is a class of flexible prior distributions for $P$. We employ the representation (2.3) with the Norm-IFPP. The advantages are as follows. We can directly estimate $M$ and $k$. This is a major difference from Miller and Harrison (2018), which estimates $M$ indirectly through the number of clusters $k$. Furthermore, an efficient Gibbs sampler can be constructed by incorporating the data augmentation with a latent Gamma random variable. This data augmentation enables us to overcome the lack of conjugacy in the categorical distribution. Instead of using the probability density function of $P_{\pi}$, we can use the density function of $h$. This is the key to building an efficient MCMC algorithm for the proposed model.

Our study is in the spirit of nonparametric Bayes in that it accounts for uncertainty in the prior distribution of the parameters by means of a random measure based on the Norm-IFPP. The details of Norm-IFPP and and the independent finite point process (IFPP) are given in Argiento and De Iorio (2022).

We illustrate the data augmentation and conditional posterior distribution of $M_{na}$. We introduce the data augmentation in models (2.1) and (2.2). This technique is employed in James et al. (2009) and Argiento and De Iorio (2022). Let $\mathcal{M}_{a}$ and $\mathcal{M}_{na}$ be the allocated and unallocated index sets, respectively. The conditional joint distribution of $\xi=(\xi_{1},\dots,\xi_{M})^{\top}$ and $S=(S_{1},\dots,S_{M})^{\top}$ given $M$ and a label vector $c$ is

where $n_{m}=\#\{i:c_{i}=m\}$. Since this equation consists of categorical distributions, conjugacy with such a distribution is required as $P_{\pi}$ for Gibbs sampling. This restriction is relaxed by using latent a random variable $U_{n}|T\sim\mathrm{Gamma}(n,T)$, where $T=\sum_{m=1}^{M}S_{m}$. In fact,

Thus, the conditional distribution of each $S_{m}$ is

If it is easy to generate random variables from the distributions in (2.4) and (2.5), an efficient Gibbs sampling algorithm can be constructed. For example, when the prior distribution of the unnormalized weight is an inverse Gaussian distribution, (2.4) and (2.5) are generalized inverse Gaussian distributions. Introducing $U_{n}$, the update of the variable $\pi$ is replaced by the update of the variable $S_{m}$. Thus, the selection of $h$ is essential in MCMC updates.

Under this data augmentation, the conditional distribution of $M_{na}$ is established in Theorem 5.1 in Argiento and De Iorio (2022). This theorem states that if $\mathcal{P}(q_{M},h,p_{0})$ follows a Norm-IFPP, then the posterior distribution of the random measure $P$, given $U_{n}$, is a superposition of a finite point process with fixed points and an IFPP. The IFPP characterizes the process of unallocated jumps, where the discrete probability distribution that serves as its parameter corresponds to the distribution of $M_{na}$. The conditional distribution of $M_{na}$ is given by

where $\theta=(\theta_{1},\dots,\theta_{n})^{\top}$, $k$ is the number of cluster (unique values of $\theta$) and $\psi$ is a Laplace transform of $h$. We sample $M_{na}$ from (2.6) in the MCMC algorithm. Then, we straightforwardly obtain $M$ by adding $k$ to $M_{na}$.

We propose a mixture of finite mixtures with normalized inverse Gaussian weights (denoted by MFM-Inv-Ga), where the notation explicitly specifies $h$ because it is essential for computation. Our proposal model only requires (2.2) to be

where $\text{Inv-Ga}(\alpha,\beta)$ is the inverse Gaussian distribution with shape parameter $\alpha>0$ and scale parameter $\beta>0$, and $\text{Norm-Inv-Ga}(\alpha_{1},\dots,\alpha_{d})$ is the normalized inverse Gaussian distribution with parameters $\alpha_{i}>0$ for $i=1,\dots,d$. The probability density function of the inverse Gaussian distribution is given by

where $\alpha>0$ is the shape parameter and $\beta>0$ is the scale parameter. The mean and variance of $S$ are given by

From Lijoi et al. (2005), the probability density function of the normalized inverse Gaussian distribution with parameters $\alpha_{i}>0,\ i=1,\dots,d$ is given by

where $\mathcal{A}_{d}(\pi_{1},\dots,\pi_{d-1})=\sum_{i=1}^{d-1}\alpha_{i}^{2}/\pi_{i}+\alpha_{d}^{2}\left(1-\sum_{i=1}^{d-1}\pi_{i}\right)^{-1}$ and $K_{m}$ is the modified Bessel functions of the third kind of order $m$ (see, also Ghosal and van der Vaart, 2017). Due to the breakdown of conjugacy with the categorical distribution, constructing an efficient Gibbs sampler is challenging. Because of the complexity of (3.3), the calculation in Miller and Harrison (2018) is intractable: constructing an MCMC algorithm based on the restaurant process by marginalizing out $M$ is extremely challenging. To overcome this difficulty, we use the data augmentation in Subsection 2.3 and the method proposed in Argiento and De Iorio (2022). As a result, it is sufficient to use not (3.3) but (3.1) and (3.8) to construct our MCMC algorithm.

We present a fast and efficient posterior computation algorithm for the proposed method. To this end, we adopt the blocked Gibbs sampling scheme in Argiento and De Iorio (2022). Let $\eta$ be a hyper-parameter of the prior distribution $q_{M}$, and $p(\cdot)$ be a joint prior density function except for $h$ and $q_{M}$. $n_{m}=\#\{i;c_{i}=m\}$ is the size of the $m$th cluster. $\psi(u)$ is the Laplace transform of $h$, which is defined by $\psi(u):=\mathbb{E}[e^{-uS_{m}}]$ for $u\geq 0$.

We summarize the algorithm in Algorithm 1. The point of this algorithm is the data augmentation through the latent variable such as $U_{n}|T\sim\mathrm{Gamma}(n,T)$, where $n$ is the sample size and $T=\sum_{m=1}^{M}S_{m}$. It is important to sample the number of empty components $M_{na}$ from (2.6) in Step 4. This step allows for direct sampling of $M$ by adding $k$ and label variables to be updated as in the finite mixture model with given $M$ in step 2. The update is the same as the telescoping sampling proposed by Frühwirth-Schnatter et al. (2021), and the method is more efficient than the classical restaurant process. However, in implementation, it is important to determine $q_{M}$ so that the series $\Psi(u,k)$ in Steps 3 and 4 can be written analytically and random variables can be easily generated from the full conditional distributions of $\eta$ and $M_{na}$. In Steps 5 and 7, if $h$ is the $\text{Inv-Ga}(\alpha,1)$, each generalized inverse Gaussian distribution (denoted by GIG) is immediately derived from (2.4) and (2.5). It is easy to generate random numbers from the GIG distribution. In Steps 6 and 7, the assigned unnormalized weights and kernel parameters are updated. When $M_{na}$ sampled in Step 4 is greater than or equal to $1$, Steps 8 and 9 are executed. Thus, when we get many empty components, it takes longer to compute.

Assume that $M-1$ follows a discrete probability distribution with the support $\{0,1,2,\dots\}$. In this paper, we consider $\text{Poisson}(\Lambda)$ ($\Lambda>0$) for $M-1$, because the constraints in Steps 3 and 4 are satisfied. In fact, from Argiento and De Iorio (2022), we have

Assuming $\Lambda\sim\text{Gamma}(a_{\Lambda},b_{\Lambda})$ for $a_{\Lambda},b_{\Lambda}>0$, the full conditional distributions of $\Lambda$ and $M_{na}$ are given by

respectively, where $\text{Shifted-Poisson}(\Lambda,t)$ is parallel shifted of the distribution $\text{Poisson}(\Lambda)$ by $t$. If we consider the negative binomial distribution as a prior distribution for $M-1$, the full conditional distributions are also obtained but learning hyper-parameters become slightly troublesome.

The choice of kernel is important, and the appropriate kernel must be selected for the purpose. In this paper, although we do not discuss the details of the selection of kernels, we present some famous kernels for the sake of completeness.

From the point of view of the interpretability of the model, the generalization performance, and the computational efficiency, it is reasonable that the number $M_{na}$ should be small. For the full conditional distribution of $M_{na}$, the following inequality holds, where $M-1|\Lambda\sim\mathrm{Poisson}(\Lambda)\ \text{and}\ \Lambda\sim\mathrm{Gamma}(a_{\Lambda},b_{\Lambda})\ \text{for}\ a_{\Lambda},b_{\Lambda}>0$.

Proposition 3.1 shows that $\psi(u)$ plays an important role in the generation of empty components. The critical difference between the inverse Gaussian and gamma distributions in estimating $M$ in the MFM is the Laplace transform. Let $\psi_{\text{Inv-Ga}}(u)$ and $\psi_{\text{Ga}}(u)$ be the Laplace transforms of $\text{Inv-Ga}(\alpha,\beta)$ and $\text{Gamma}(\gamma,\eta)$, respectively. The we have

Laplace transforms $\psi_{\text{Inv-Ga}}$ and $\psi_{\text{Ga}}$ are decreasing functions with respect to $u$. The former has exponential decay, while the latter is polynomial. Figure 1 shows the graphs of the Laplace transforms of inverse Gaussian and gamma distributions when the shape parameters are $\alpha,\gamma=1.0,0.2,10^{-1},10^{-2},10^{-3}$ and the scale parameters are $1$. It can be seen that $\psi_{\text{Inv-Ga}}$ decreases much faster than $\psi_{\text{Ga}}$. The speed of this decrease has a significant impact on the appearance of the empty component in estimating the number of components.

Inequalities (3.6) and (3.7) are not for marginal posterior distributions, but for full conditional distributions of $M_{na}$. This is due to the fact that the marginal posterior distributions of $U_{n}$ are analytically intractable.

In this subsection, we illustrate the performance of clustering and inference for the number of components using artificial and real data.