Penalized Maximum Likelihood Estimator for Mixture of von Mises-Fisher Distributions

In many areas of statistical modelling, data are represented as directions or unit length vectors (Mardia, 1972; Jupp, 1995; Mardia and Jupp, 2000). The analysis of directional data has attracted much research interests in various disciplines, from hydrology (Chen et al., 2013) to biology (Boomsma et al., 2006), and from image analysis (Zhe et al., 2019) to text mining (Banerjee et al., 2005). The von Mises-Fisher (vMF) distribution is one of the most commonly used distributions to model data distributed on the surface of the unit hypersphere (Fisher, 1953; Mardia and Jupp, 2000). The vMF distribution has been applied successfully in many domains (e.g., (Sinkkonen and Kaski, 2002; Mcgraw et al., 2006; Bangert et al., 2010)).

A mixture of vMF distributions (Banerjee et al., 2005) assumes that each observation is drawn from one of the $p$ vMF distributions. Applications of the vMF mixture model are diverse, including image analysis (Mcgraw et al., 2006) and text mining (Banerjee et al., 2005). More recently, it has been shown that the vMF mixture model can approximate any continuous density function on the unit hypersphere to arbitrary degrees of accuracy given sufficient numbers of mixture component (Ng and Kwong, 2021).

Various estimation strategies have been developed to perform model estimation, including the maximum likelihood approach (Banerjee et al., 2005) and Bayesian methods (Bagchi and Kadane, 1991; Taghia et al., 2014). The maximum likelihood approach, which is typically performed using the Expectation-Maximization (EM) algorithm (Dempster et al., 1977; Banerjee et al., 2005), is among the most popular approach to parameter estimation. However, as we show in Section 3 the likelihood function is unbounded from above, and consequently a global maximum likelihood estimate (MLE) fails to exist.

The unboundedness of likelihood function occurs in various mixture modelling context, particularly for mixture models with location-scale family distributions including the mixture of normal distributions (Ciuperca et al., 2003; Chen et al., 2008) and the mixture of Gamma distributions (Chen et al., 2016). Various approaches have been developed in order to tackle the likelihood degeneracy problem, including resorting to local estimates (Peters et al., 1978), compactification of the parameter space (Redner, 1981), and constrained maximization of the likelihood function (Hathaway, 1985).

An alternative solution to the likelihood degeneracy problem is to penalized the likelihood function such that the resulting penalized likelihood function is bounded and the existence of the penalized MLE is guaranteed. The approach of penalized maximum likelihood was applied to normal mixture models (Ciuperca et al., 2003; Chen et al., 2008), two-parameter Gamma mixture models (Chen et al., 2016). A penalized likelihood approach has also a Bayesian interpretation (Goodd and Gaskins, 1971; Ciuperca et al., 2003), whereby the penalized likelihood function corresponds to a posterior density and the penalized maximum likelihood solution to the maximum a posterior estimate.

Previously, the penalized maximum likelihood approach is applied to study the mixture of von-Mises distributions (Chen et al., 2007) where consistency results were obtained. The von-Mises distribution is a special case of the von Mises-Fisher distribution defined on the circle. We generalize the results in Chen et al. (2007) to the arbitrary dimensional sphere. The consistency proof in Chen et al. (2007) relies heavily on the univariate properties of the von Mises distribution and generalization of the result to higher dimensions is not straightforward. In this paper we prove a few useful technical lemmas before proving the main results. To handle the non-identifiability of the mixture models, we use the framework of Redner (1981) to obtain consistency in the quotient space.

In this paper, we also consider the penalized likelihood approach to tackle the problem of likelihood unboundedness for the mixture of vMF distributions. We incorporate a penalty term into the likelihood function and maximize the resulting penalized likelihood function. We study conditions on the penalty function to ensure consistency of the penalized maximum likelihood estimator (PMLE). We develop an Expectation-Maximization algorithm to perform model estimation based on the penalized likelihood function. The rest of the paper is structured as follows. Section 2 introduces the background on vMF mixtures and key notations used in the subsequent sections. The problem of likelihood degeneracy is formally presented in Section 3. Section 4 develops the penalized maximum likelihood approach and discusses conditions on the penalty function to ensure strong consistency of the resulting estimator. An Expectation-Maximization algorithm is also developed in Section 4, and its performance is examined in Section 5. Section 6 illustrate the proposed EM algorithm using a data application. We conclude the paper with a discussion section.

The probability density function of a $d$-dimensional vMF distribution is given by:

where $x\in\mathbb{S}^{d-1}$ is a $d$-dimensional unit vector (i.e. $\mathbb{S}^{d-1}:=\{\mathbf{x}\in\mathbb{R}^{d}:||\mathbf{x}||=1\}$ where $||\cdot||$ is the $L_{2}$ norm), $\boldsymbol{\mu}\in\mathbb{S}^{d-1}$ is the mean direction, and $\kappa\geq 0$ is the concentration parameter. The normalizing constant $c_{d}(\kappa)$ has the form

where $I_{r}(\cdot)$ is the modified Bessel function of the first kind of order $r$. The vMF distribution becomes increasingly concentrated at the mean direction $\mu$ as the concentration parameter $\kappa$ increases. The case $\kappa=0$ corresponds to the uniform distribution on $\mathbb{S}^{d-1}$.

The probability density function of a mixture of vMF distributions with $p$ mixture components can be expressed as

where $(\pi_{1},\ldots,\pi_{p})$ is the mixing proportions, $(\boldsymbol{\mu}_{k},\kappa_{k})$ are the parameters for the $k$th component of the mixture, and $f(\cdot;\boldsymbol{\mu}_{k},\kappa_{k})$ is the vMF density function defined in (1).

We let $\Theta:=\{\theta\equiv(\boldsymbol{\mu},\kappa):\boldsymbol{\mu}\in\mathbb{S}^{d-1},\kappa\geq 0\}$ be the parameter space of the vMF distribution, with the metric $\rho(\cdot,\cdot)$ defined as

for $\theta_{1}=(\boldsymbol{\mu}_{1},\kappa_{1}),\theta_{2}=(\boldsymbol{\mu}_{2},\kappa_{2})$. For any $\theta=(\boldsymbol{\mu},\kappa)\in\Theta$, we write $f_{\theta}(\cdot):=f(\cdot;\boldsymbol{\mu},\kappa)$ for the density function and $\gamma_{\theta}$ for the corresponding measure. The space of mixing probabilities is denoted by $\Pi:=\{(\pi_{1},\ldots,\pi_{p}):\sum_{i=1}^{p}\pi_{p}=1,\pi_{k}\geq 0,k=1,\ldots,p\}$. A $p$-component mixture of vMF distributions can be expressed as $\gamma=\sum_{k=1}^{p}\pi_{k}\gamma_{\theta_{k}}$ where $(\pi_{1},\ldots,\pi_{p})\in\Pi$ and $(\theta_{1},\ldots,\theta_{p})\in\Theta^{p}$, and where $\Theta^{p}=\Theta\times\cdots\times\Theta$ is the product of the parameter spaces. We define the product space $\Gamma:=\Pi\times\Theta^{p}$, and we slightly abuse notations to let $\gamma$ denote both the mixing measure and the parameters in $\Gamma$. While $\Gamma$ is a natural parameterization of the family of mixture of vMF distributions, elements of $\Gamma$ are not identifiable. Thus, we let $\tilde{\Gamma}$ be the quotient topological space obtained from $\Gamma$ by identifying all parameters $(\pi_{1},\ldots,\pi_{p},\theta_{1},\ldots,\theta_{p})$ such that their corresponding densities are equal (almost) everywhere. For the rest of the paper, we assume that the number of mixture components $p$ is known.

We investigate the likelihood degeneracy problem of the vMF mixture model in this section. For any observations generated from a vMF mixture model with two or more mixture components, we show that the resulting likelihood function on the parameter space $\Gamma$ is unbounded above. As discussed in Section 1, likelihood degeneracy is a common problem for mixture models with location-scale distributions, including the normal mixtures. In the case of normal mixture distributions, one can show that by letting the mean parameter of a mixture component equal to one of the observations and letting the variance of the same mixture component converge to zero while holding other parameters fixed, the likelihood function diverges to positive infinity (Chen et al., 2008).

For the vMF mixture distributions, the likelihood unboundedness can be best understood in the special case of $x\in\mathbb{S}^{1}$, or the mixture of von Mises distributions. The von Mises distribution, also known as the circular normal distribution, approaches a normal distribution with large concentration parameter $\kappa$:

with $\sigma^{2}=1/\kappa$, and the approximation converges uniformly as $\kappa$ goes to infinity. Therefore, the likelihood function of a mixture of von Mises distributions diverges to infinity by letting the mean parameter of a mixture component equal to one of the observations and letting the concentration parameter diverges to infinity.

We now consider the general case of the vMF mixture models. Let ${\cal X}=\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\}$ be the observations generated from a mixture of vMF distributions with density function $\sum_{k=1}^{p}\pi_{k}f_{\theta_{k}}(\cdot)$ where $\theta_{k}=(\boldsymbol{\mu}_{k},\kappa_{k})$. The likelihood function can be expressed as:

where $\boldsymbol{\theta}=(\theta_{1},\ldots,\theta_{p})=((\boldsymbol{\mu}_{1},\kappa_{1}),\ldots,(\boldsymbol{\mu}_{p},\kappa_{p}))$ and $\boldsymbol{\pi}=(\pi_{1},\ldots,\pi_{p})$. We can show that by letting the mean direction $\boldsymbol{\mu}_{k}$ of one of the mixture components equals to an arbitrary observation and letting the corresponding concentration parameter $\kappa_{k}$ goes to infinity, the resulting likelihood function diverges.

The proof of Theorem 3.1 is provided in the Appendix. The unboundedness of the likelihood function on the parameter space implies that the maximum likelihood estimator is not well defined.

We perform simulation studies to investigate the performance of the proposed EM algorithm for maximizing the penalized likelihood function. We generate data from the mixture of vMF distributions with two and three mixture components and with dimensions $d=2,3,4$. For each model, data are generated with increasing samples sizes to assess the convergence of the estimated parameters toward the true parameters. The concentration parameters $\boldsymbol{\kappa}$ and the mixing proportions $\boldsymbol{\pi}$ are pre-specficied whereas the mean directions $\boldsymbol{\mu}$ are drawn from the uniform distribution on the surface of the unit hypersphere.

For the two mixture components model, we specify the mixing proportions as $\boldsymbol{\pi}=(0.5,0.5)$ and the concentration parameters $\boldsymbol{\kappa}=(10,1)$. For the model with three mixture components, we set $\boldsymbol{\pi}=(0.4,0.3,0.3)$ and $\boldsymbol{\kappa}=(10,5,1)$. For illustrative purpose, we consider the penalty function $\tilde{p}_{n}(\kappa_{l})=-(1/n)\kappa_{l}$. For each combination of dimension $d$ and sample size $n$, we simulate 500 random samples from the model and the EM algorithm developed in Section 4.3 is used to obtain the parameter estimates. We measure the distance between the estimated parameters and the true parameters for each random sample. For the mean direction parameters $\boldsymbol{\mu}$, the distance is measured using the metric $d(\mathbf{x},\mathbf{y})=\mbox{arccos}(\mathbf{x}^{T}\mathbf{y})$.

Simulation results for the two and three mixture cases are presented in Table 1 and 2, respectively. The average distance and the standard deviation between the true and the estimated parameters from 500 replications are reported. We observe that the estimated parameters converge to the true parameter as $n$ increases. We notice that the mean direction parameter can be estimated with higher precision when the corresponding concentration parameter is large. This is expected since observations are more closely clustered with a large concentration parameter.

Table 3 and 4 show the number of degeneracies when running the EM algorithm for computing the ordinary MLE for mixture of vMF distributions. Observations are generated from mixture of vMF distributions with one mixture component for Table 3 and with two mixture components for Table 4. We vary the dimension of the data from $d=3$ to $d=4$ and the sample size from $n=100$ to $n=500$. Mixtures of vMF distributions with $p=2,3,4,5$ components are fitted to the generated data. We compute the ordinary MLE using the EM algorithm and record the number of times that the EM fails to converge from 1000 simulation runs. The EM algorithm is considered fail to converge if one of the concentration parameters becomes exceedingly large (greater than $10^{10}$). From Table 3 and 4, the EM algorithm tends to fail to converge with smaller sample sizes. We also note that when the fitted model has a larger number of mixture components $p$, the EM algorithm is more likely to fail to converge.

We illustrate the EM algorithm for maximum penalized log-likelihood using the household data set from R package HSAUR3. The data set contains the household expenditures of 20 single men and 20 single women on four commodity group. As in Hornik and Grün (2014), we will also focus on the three commodity groups (housing, food and service). The EM algorithms for ordinary MLE and for the penalized MLE with 2 and 3 mixture components are fitted to the data. The results are shown in Table 5 and 6, respectively, where the estimated parameters $\hat{\boldsymbol{\pi}},\hat{\boldsymbol{\mu}},\hat{\boldsymbol{\kappa}}$ are shown for all cases. The estimated paramters for the MLE and for the penalized MLE are very similar for both $p=2$ and $p=3$. The log-likelihood evaluated at the MLE is slightly larger than the penalized log-likelihood evaluated at the penalized MLE. More interestingly, we observe that for each case the largest concentration parameter obtained under the penalized MLE is smaller than that obtained under the MLE. This behavior suggests that the incorporation of a penalty function pulls back the estimate of largest concentration parameter towards 0 and prevents the divergence of the likelihood function.

In this paper we considered a penalized maximum likelihood approach to the estimation of the mixture of vMF distributions. By incorporating a suitable penalty function, we showed that the resulting penalized MLE is strongly consistent. An EM algorithm was derived to maximize the penalized likelihood function, and its performance and behavior were examined using simulation studies and a data application. The techniques used in this work to prove consistency could be applicable to study other mixture models for spherical observations.

On behalf of all authors, the corresponding author states that there is no conflict of interest.