On Model-Based Clustering of Directional Data with Heavy Tails

Abstract

Directional statistics deals with data that can be naturally expressed in the form of vector directions. The von Mises-Fisher distribution is one of the most fundamental parametric models to describe directional data. Mixtures of von Mises-Fisher distributions represent a popular approach to handling heterogeneous populations. However, components of such models can be affected by the presence of mild outliers or cluster tails heavier than what can be accommodated by means of a von Mises-Fisher distribution. To relax these model limitations, a mixture of contaminated von Mises-Fisher distributions is proposed. The performance of the proposed methodology is tested on synthetic data and applied to text and genetics data. The obtained results demonstrate the importance of the proposed procedure and its superiority over the traditional mixture of von Mises-Fisher distributions in the presence of heavy tails.

Appendix

1.1 M-Step Derivation

Based on the Q-function provided in (3), closed-form expressions for parameters $\pi _g$, $\delta _g$, and $\varvec{\mu }_g$ can be obtained. Parameters $\kappa _g$ and $\lambda _g$ need to be estimated numerically. For parameter $\kappa _g$, we obtain

$$\begin{aligned} \frac{\partial Q(\kappa _g)}{\partial \kappa _g} = \sum _{i=1}^{n} \tau ^{(s)}_{ig} \bigg [v_{i|g}^{(s)} \frac{c_p^{\prime } (\kappa _g)}{c_p(\kappa _g)} + v^{(s)}_{i|g} \varvec{\mu }_g^{(s)\top } \varvec{x}_i + (1-v^{(s)}_{i|g}) \frac{c_p^\prime (\lambda _g^{(s-1)} \kappa _g) }{c_p(\lambda ^{(s-1)}_g \kappa _g)} + (1-v^{(s)}_{i|g}) \lambda _g^{(s-1)} \varvec{\mu }_g^{(s) \top } \varvec{x}_i\bigg ]. \end{aligned}$$

(10)

As $c_p(\kappa _g) = (2\pi )^{-\frac{p}{2}} \kappa _g^{\frac{p}{2} - 1} I^{-1}_{\frac{p}{2}-1}(\kappa _g)$, by denoting $r=\frac{p}{2}-1$ we obtain $c_p(\kappa _g) = \frac{\kappa _g^r}{(2\pi )^{r+1} I_r(\kappa _g)}$ and $c_p^{\prime } (\kappa _g) = \frac{r \kappa _g^{r-1} I_r(\kappa _g) - \kappa _g^r I_r^{\prime }(\kappa _g)}{(2\pi )^{r+1} I^2_r(\kappa _g)}$. Taking into account the fact that $I_r^{\prime } (\kappa _g) = \frac{r}{\kappa _g} I_r (\kappa _g) + I_{r+1} (\kappa _g)$, it yields

$$\begin{aligned} \frac{c_p^{\prime } (\kappa _g)}{c_p(\kappa _g)} = \frac{\big [r\kappa _g^{r-1} I_r(\kappa _g) - \kappa _g^r (\frac{r}{\kappa _g} I_r(\kappa _g) + I_{r+1}(\kappa _g))\big ] (2\pi )^{r+1} I_r(\kappa _g)}{(2\pi )^{r+1} I^2_r(\kappa _g) \kappa _g^r} = -\frac{I_{r+1}(\kappa _g)}{I_r(\kappa _g)}. \end{aligned}$$

(11)

Similarly, it follows that

$$\begin{aligned} \frac{c_p^\prime (\lambda _g^{(s-1)} \kappa _g) }{c_p(\lambda _g^{(s-1)} \kappa _g)} = -\lambda _g^{(s-1)} \frac{I_{r+1}(\lambda _g^{(s-1)} \kappa _g)}{I_r(\lambda _g^{(s-1)} \kappa _g)}. \end{aligned}$$

(12)

Employing the Newton–Raphson method for the maximization of the Q-function with respect to $\kappa _g$, we have $\kappa _g^{(s)[j]} = \kappa _g^{(s)[j-1]}-\frac{U(\kappa _g^{(s)[j-1]})}{U^\prime (\kappa _g^{(s)[j-1]})}$, where j is the iteration number of the Newton–Raphson algorithm and $U(\kappa _g)=\frac{\partial Q(\kappa _g)}{\partial \kappa _g}$ is given in (10). We proceed to obtain the derivative of $U(\kappa _g)$ as follows below:

$$\begin{aligned} U^{\prime } (\kappa _g) = \sum _{i=1}^{n} \tau _{ig}^{(s)} v_{i|g}^{(s)} \frac{\partial }{\partial \kappa _g}\bigg (-\frac{I_{r+1}(\kappa _g)}{I_r(\kappa _g)}\bigg ) -\lambda _g^{(s-1)} \sum _{i=1}^{n} \tau _{ig}^{(s)} (1-v_{i|g}^{(s)}) \frac{\partial }{\partial \kappa _g}\bigg (\frac{I_{r+1}(\lambda _g^{(s-1)} \kappa _g)}{I_r(\lambda _g^{(s-1)} \kappa _g)}\bigg ), \end{aligned}$$

(13)

where

$$\begin{aligned} \begin{aligned} \frac{\partial }{\partial \kappa _g}\bigg (\frac{I_{r+1}(\kappa _g)}{I_r(\kappa _g)}\bigg )&= \frac{[\frac{r+1}{\kappa _g} I_{r+1}(\kappa _g) + I_{r+2}(\kappa _g)] I_r(\kappa _g) - I_{r+1}(\kappa _g) [\frac{r}{\kappa _g} I_r(\kappa _g) + I_{r+1}(\kappa _g)]}{I_r^2(\kappa _g)} \\&= \frac{I_{r+2}(\kappa _g)}{I_r(\kappa _g)} + \frac{1}{\kappa _g} \frac{I_{r+1}(\kappa _g)}{I_r(\kappa _g)} - \frac{I_{r+1}^2(\kappa _g)}{I_r^2(\kappa _g)} \end{aligned} \end{aligned}$$

(14)

and

$$\begin{aligned} \begin{aligned} \frac{\partial }{\partial \kappa _g}&\bigg (\frac{I_{r+1}(\lambda _g^{(s-1)} \kappa _g)}{I_r(\lambda _g^{(s-1)} \kappa _g)}\bigg ) = \frac{ \frac{\partial I_{r+1}(\lambda _g^{(s-1)} \kappa _g)}{\partial \kappa _g} I_r(\lambda _g^{(s-1)} \kappa _g) - I_{r+1}(\lambda _g^{(s-1)} \kappa _g) \frac{\partial I_r(\lambda _g^{(s-1)} \kappa _g)}{\partial \kappa _g}}{I_r^2(\lambda _g^{(s-1)} \kappa _g)} \\&= \frac{\big [\frac{r+1}{\lambda _g^{(s-1)} \kappa _g} I_{r+1}(\lambda _g^{(s-1)} \kappa _g) + I_{r+2}(\lambda _g^{(s-1)} \kappa _g)\big ] \lambda _g^{(s-1)} I_r(\lambda _g^{(s-1)} \kappa _g)}{I_r^2(\lambda _g^{(s-1)} \kappa _g)} \\&\quad - \frac{I_{r+1}(\lambda _g^{(s-1)} \kappa _g) \big [\frac{r}{\lambda _g^{(s-1)} \kappa _g} I_r(\lambda _g^{(s-1)} \kappa _g) + I_{r+1}(\lambda _g^{(s-1)} \kappa _g)\big ] \lambda _g^{(s-1)}}{I_r^2(\lambda _g^{(s-1)} \kappa _g)}\\&= \lambda _g^{(s-1)} \frac{I_{r+2}(\lambda _g^{(s-1)} \kappa _g)}{I_r(\lambda _g^{(s-1)} \kappa _g)} + \frac{1}{\kappa _g} \frac{I_{r+1}(\lambda _g^{(s-1)} \kappa _g)}{I_r(\lambda _g^{(s-1)} \kappa _g)}- \lambda _g^{(s-1)} \frac{I_{r+1}^2(\lambda _g^{(s-1)} \kappa _g)}{I_r^2(\lambda _g^{(s-1)} \kappa _g)}. \end{aligned} \end{aligned}$$

(15)

Substituting (11)-(15) into (10) yields the result presented in (5) and (6).

Similarly, for $\lambda _g$ we have

$$\begin{aligned} \frac{\partial Q(\lambda _g)}{\partial \lambda _g} = \sum _{i=1}^{n} \tau _{ig}^{(s)} (1-v_{i|g}^{(s)})\bigg [\frac{1}{c_p(\lambda _g \kappa _g^{(s)})}\frac{\partial c_p(\lambda _g \kappa _g^{(s)}) }{\partial \lambda _g } + \kappa _g^{(s)} \varvec{\mu }_g^{(s) \top } \varvec{x}_i\bigg ]. \end{aligned}$$

Now, noticing that $\frac{1}{c_p(\lambda _g \kappa _g^{(s)})}\frac{\partial c_p(\lambda _g \kappa _g^{(s)}) }{\partial \lambda _g } = -\kappa _g^{(s)} \frac{I_{r+1}(\lambda _g \kappa _g^{(s)})}{I_r (\lambda _g \kappa _g^{(s)})}$ and making use of

$$\begin{aligned} \frac{\partial }{\partial \lambda _g}\bigg (\frac{I_{r+1}(\lambda _g \kappa _g^{(s)})}{I_r(\lambda _g \kappa _g^{(s)})}\bigg ) = \kappa _g^{(s)} \frac{I_{r+2}(\lambda _g \kappa _g^{(s)})}{I_r(\lambda _g \kappa _g^{(s)})} + \frac{1}{\lambda _g} \frac{I_{r+1}(\lambda _g \kappa _g^{(s)})}{I_r(\lambda _g \kappa _g^{(s)})}- \kappa _g^{(s)} \frac{I_{r+1}^2(\lambda _g \kappa _g^{(s)})}{I_r^2(\lambda _g \kappa _g^{(s)})}, \end{aligned}$$

we obtain (7) and (8) for $V(\lambda _g)=\frac{\partial Q(\lambda _g)}{\partial \lambda _g}$.

1.2 $c_p(\kappa )$ Approximation

Let $H_p(\kappa )=\frac{1}{c_{\frac{p}{2}-1}(\kappa )}$. Then, based on the Amos-type bounds obtain

$$\begin{aligned} S_{p+\frac{1}{2},p+\frac{3}{2}} (\kappa ) \le \log (H_{p}(\kappa )) \le \min (S_{p,p+2}(\kappa ), S_{p+\frac{1}{2},\sqrt{(p+\frac{1}{2})(p+\frac{3}{2})}}(\kappa )), \end{aligned}$$

where $S_{\alpha ,\beta }(\kappa )=\sqrt{\kappa ^2+\beta ^2}-\alpha \log (\alpha +\sqrt{\kappa ^2+\beta ^2})-\beta +\alpha \log (\alpha +\beta )$. Let

$$\begin{aligned} L_{p}(\kappa )=S_{p+\frac{1}{2}, \sqrt{(p+\frac{1}{2})(p+\frac{3}{2})}}(\kappa ) + S_{p,p+2}(\min (\kappa ,\kappa _{p})-S_{p+\frac{1}{2}, \sqrt{(p+\frac{1}{2})(p+\frac{3}{2})}}(\min (\kappa ,\kappa _{p})) \end{aligned}$$

(16)

denote an approximation for $\log (H_{p}(\kappa ))$, where $\kappa _{p}=\sqrt{(3p+\frac{11}{2})(p+\frac{3}{2})}$. First, calculate $L_{p}(\kappa )$ based on (16). For a threshold value $\theta =700$ chosen in such a way that $e^{\theta }$ does not overflow and $e^{-\theta }$ does not underflow at the first step, calculate the logarithm of $H_{p}(\kappa )$ directly if $L_{p}(\kappa ) \le \theta -\frac{1}{2}$. If $\theta - \frac{1}{2} < L_{p}(\kappa ) \le 2\theta -1$, use the approximation

$$\begin{aligned} \log H_{p}(\kappa )=\frac{L_{p}(\kappa )}{2} + \log \Big (\sum _{j=0}^{\infty } \frac{e^{\frac{-L_{p}(\kappa )}{2}} \Gamma (p+1) }{\Gamma (p+1+j)} \frac{(\frac{\kappa ^2}{4})^j}{j!}\Big ). \end{aligned}$$

Otherwise, use $L_{p}(\kappa )$ as an approximation for $\log (H_{p}(\kappa ))$. The reader is referred to Hornik and Grün (2014) for more details.