Finite mixtures of unimodal beta and gamma densities and the $k$-bumps algorithm

Abstract

This paper addresses the problem of estimating a density, with either a compact support or a support bounded at only one end, exploiting a general and natural form of a finite mixture of distributions. Due to the importance of the concept of multimodality in the mixture framework, unimodal beta and gamma densities are used as mixture components, leading to a flexible modeling approach. Accordingly, a mode-based parameterization of the components is provided. A partitional clustering method, named $k$-bumps, is also proposed; it is used as an ad hoc initialization strategy in the EM algorithm to obtain the maximum likelihood estimation of the mixture parameters. The performance of the $k$-bumps algorithm as an initialization tool, in comparison to other common initialization strategies, is evaluated through some simulation experiments. Finally, two real applications are presented.

Appendices

Parameterization genesis

If a density function $f$ is chosen to belong to the Pearson system, then it is the solution of the differential equation

$$\begin{aligned} \frac{df\left(x\right)}{dx}=-\frac{\left(x-m\right)f\left(x\right)}{c_0+c_1x+c_2x^2}. \end{aligned}$$

(18)

It is clear that some of the solutions to (18) have a single mode ($df/dx=0$ at $x=m$) and smooth contact with the horizontal axis ($df/dx=0$ when $f\left(x\right)=0$).

The shape of $f$ depends on the set of parameters $\left(m,c_0,c_1,c_2\right)$, and the form of the solution of (18) evidently depends on the nature of the roots of the equation

$$\begin{aligned} c_0+c_1x+c_2x^2=0. \end{aligned}$$

(19)

The classes of unimodal gamma and beta densities, illustrated in Sect. 2, arise from a convenient choice of these roots.

1.1 Unimodal gamma densities

Consider $c_2=0$ and $c_1>0$. Thus, Eq. (18) becomes

$$\begin{aligned} \frac{df\left(x\right)}{dx}=-\frac{\left(x-m\right)f\left(x\right)}{c_0+c_1x}=\left(-\frac{1}{c_1}+\frac{m+c_0/c_1}{c_0+c_1x}\right)f\left(x\right), \end{aligned}$$

(20)

in which

$$\begin{aligned} f\left(x\right)=C\left(x+\frac{c_0}{c_1}\right)^{\frac{1}{c_1}\left(m+\frac{c_0}{c_1}\right)}e^{-\frac{x}{c_1}},\qquad -\frac{c_0}{c_1}\le x<\infty . \end{aligned}$$

(21)

In order to make (21) a density function,

$$\begin{aligned} C=\left\{ c_1^{\frac{1}{c_1}\left(m+\frac{c_0}{c_1}\right)+1}e^{\frac{c_0}{c_1^2}}\mathrm \Gamma \left[\frac{1}{c_1}\left(m+\frac{c_0}{c_1}\right)+1\right]\right\} ^{-1}, \end{aligned}$$

so that the result is a gamma distribution (see Johnson and Kotz 1970a, Chapter 17). Equation (2) is obtained by setting $-c_0/c_1=a$ and $c_1=v$.

1.2 Unimodal beta densities

Suppose that both the roots of (19) are real. Denoting these roots as $a$ and $b$, with $a<b$, it follows that

$$\begin{aligned} c_0+c_1x+c_2x^2=-c_2\left(x-a\right)\left(b-x\right); \end{aligned}$$

consequently, Eq. (18) becomes

$$\begin{aligned} \frac{df\left(x\right)}{dx}=\frac{\left(x-m\right)f\left(x\right)}{c_2\left(x-a\right)\left(b-x\right)}=\frac{1}{c_2\left(b-a\right)}\left(\frac{a-m}{x-a}+\frac{b-m}{b-x}\right)f\left(x\right), \end{aligned}$$

(22)

in which

$$\begin{aligned} f\left(x\right)=C\left(x-a\right)^{\frac{a-m}{c_2\left(b-a\right)}}\left(b-x\right)^{\frac{m-b}{c_2\left(b-a\right)}},\qquad a\le x\le b. \end{aligned}$$

(23)

In order to make (23) a density function,

$$\begin{aligned} C=\left\{ \left(b-a\right)^{\frac{c_2-1}{c_2}}\mathrm B \left[\frac{a-m}{c_2\left(b-a\right)}+1,\frac{m-b}{c_2\left(b-a\right)}+1\right]\right\} ^{-1}, \end{aligned}$$

so that the result is a beta distribution (see Johnson and Kotz 1970b, Chapter 24). If a unimodal beta density must be considered, it is necessary that $c_2\le 0$. Equation (6) is obtained by setting $-c_2=v$ in (23).

Details on the EM algorithm

Here we attempt to make explicit the derivatives in (14) for both gamma and beta densities parameterized according to (2) and (6), respectively. We recall that the resulting ML-estimates do not have a closed-form expression and can only be computed numerically, with the aid of an iterative algorithm; such numerical methods are available in most computer software, such as Mathematica and R.

In detail, for the gamma density in (2) we have

$$\begin{aligned} \displaystyle \frac{\partial \ln f \left(x_i;m_j,v_j\right)}{\partial m_j} = \displaystyle \frac{1}{v_j}\left[\ln \left(x_i-a\right)-\ln v_j-\psi \left(\frac{m_j-a}{v_j}+1\right)\right] \end{aligned}$$

and

$$\begin{aligned} \frac{\partial \ln f \left(x_i;m_j,v_j\right)}{\partial v_j}&= \frac{1}{v_j^2}\left\{ \left(m_j-a\right)\left[\ln v_j+\psi \left(\displaystyle \frac{m_j-a}{v_j}+1\right)-\ln \left(x_i-a\right)\right]+\right.\\&-\left(m_j+v_j\right)+x_i\biggr \}, \end{aligned}$$

where $\psi \left(\cdot \right)$ is the digamma function. In the same way, for the beta density in (6) we have

$$\begin{aligned} \frac{\partial \ln f \left(x_i;m_j,v_j\right)}{\partial m_j}&= \displaystyle \frac{1}{v_j\left(b-a\right)}\left\{ \left[\psi \left(\frac{b-m_j}{v_j\left(b-a\right)}+1\right)-\psi \left(\frac{m_j-a}{v_j\left(b-a\right)}+1\right)\right]+\right.\\&\quad +\ln \left(x_i-a\right)-\ln \left(b-x_i\right)\biggr \}, \end{aligned}$$

and

$$\begin{aligned} \frac{\partial \ln f\left(x_i;m_j,v_j\right)}{\partial v_j}&= \frac{1}{v_j^2\left(b-a\right)}\left\{ \left(b-a\right)\left[\ln \left(b-a\right)-\psi \left(\displaystyle \frac{2v_j+1}{v_j}\right)\right]+\right.\\&\quad +\left[\left(m_j-a\right)\psi \quad \left(\frac{m_j-a}{v_j\left(b-a\right)}+1\right)+\left(b-m_j\right)\psi \left(\displaystyle \frac{b-m_j}{v_j\left(b-a\right)}+1\right)\right]+\\&\quad -\left(m_j-a\right)\ln \left(x_i-a\right)-\left(b-m_j\right)\ln \left(b-x_i\right)\biggr \}. \end{aligned}$$