Robust Fitting of Mixture Models using Weighted Complete Estimating Equations

Mixture modeling (McLachlan and Peel, 2000) is very popular for distribution estimation, regression, and model-based clustering by taking account of potential subgroup structures of data. Typically, such mixture models are fitted by the maximum likelihood method using the well-known Expectation-Maximization (EM) algorithm. However, data often contain outliers and the maximum likelihood method can be highly affected by them. The presence of outliers would result in the biased and inefficient statistical inference of the parameters of interest and would thus make recovering the underlying clustering structure of the data very difficult. A typical remedy for this problem is to replace the normal component distributions with heavy-tailed components. For example, Bagnato et al. (2017), Forbes and Wraith (2014), Peel and McLachlan (2004), Punzo and Tortora (2021), Sun et al. (2010), and Zhang and Liang (2010) used symmetric component distributions, Basso et al. (2010), Lee and McLachlan (2014), Lin (2010), Morris et al. (2019), and Wang et al. (2009) used skewed distributions, and Galimberti and Soffritti (2014), Ingrassia et al. (2014), Mazza and Punzo (2020), Punzo et al. (2017), Song et al. (2014), Yao et al. (2014), and Zarei et al. (2019) considered robustifying the mixture of regressions. However, this approach cannot distinguish non-outliers from outliers straightforwardly because it fits a mixture model to all the observations, including outliers. Apart from using heavy-tailed distributions, a class of generalized likelihood is an effective alternative tool for the robust fitting of statistical models. In the context of mixture modeling, Fujisawa and Eguchi (2006) employed density power divergence (Basu et al., 1998) for robust estimation of Gaussian mixture models. However, a direct application of the divergence may suffer from computational problems because the objective functions may contain intractable integrals. For example, the objective function under the density power divergence includes the integral $\int f(x;\theta)^{\gamma}dx$, where $f(x;\theta)$ is the density function with the parameter $\theta$ and $\gamma$ is the tuning parameter. When $f(x;\theta)$ is the density function of a mixture model, the integral becomes intractable. Therefore, these approaches are not necessarily appealing in practice, despite their desirable robustness properties.

In this paper, we develop a new approach to robust mixture modeling by newly introducing the idea of the weighted complete estimating equations (WCE). Complete estimating equations appear in the EM algorithm and are conditioned on the latent membership variables of the cluster assignment. Our WCE introduces weights to the complete estimating equations such that the contribution of outliers to the complete estimating equations is automatically downweighted. The weight is defined based on the assumed component distributions. By conditioning the membership variables in WCE, only the weighted estimating equations for each component distribution must be considered separately instead of directly using the weighted estimating equations for a mixture model. Since the derived WCE depends on the unknown latent membership variables, these latent variables are augmented via their posterior expectations to solve the WCE, which provides an expectation-estimating-equation (EEE) algorithm (Elashoff and Ryan, 2004). The proposed EEE algorithm is general and can be applied to various mixture models. The proposed WCE method is then applied to three types of mixture models: a multivariate Gaussian mixture, a mixture of experts (Jacobs et al., 1991), and a multivariate skew-normal mixture (Lin et al., 2007; Lin, 2009). For the multivariate Gaussian mixture, the updating steps of the proposed EEE algorithm are obtained in closed form without requiring either numerical integration or optimization steps. A similar algorithm can be derived for a mixture of experts when the component distributions are Gaussian. Moreover, for the multivariate skew normal mixture, by introducing additional latent variables based on the stochastic representation of the skew normal distribution, a novel EEE algorithm can be obtained as a slight extension of the proposed general EEE algorithm in which all the updating steps proceed analytically.

The proposed framework of the weighted complete estimating equations unifies several existing approaches regarding the choice of the weight function, such as hard trimming (e.g. Garcia-Escudero et al., 2008), soft trimming (e.g. Campbell, 1984; Farcomeni and Greco, 2015; Greco and Agostinelli, 2020) and their hybrid (e.g. Farcomeni and Punzo, 2020). Among several designs for the weight function, this study adopts the density weight, which leads to theoretically valid and computationally tractable robust estimation procedures not only for the Gaussian mixture but also for a variety of mixture models.

As a related work, Greco and Agostinelli (2020) employed a similar idea using the weighted likelihood for the robust fitting of the multivariate Gaussian mixture. They constructed weights by first obtaining a kernel density estimate for the Mahalanobis distance concerning the component parameters and then computing the Pearson residuals. Since their method requires knowledge of the distribution of the quadratic form of the residuals, their approach cannot be directly extended to other mixture models. On the other hand, the proposed method can be applied to various mixture models, as in the proposed weight design, and the weights are automatically determined by the component-wise density. It should be noted that introducing the density weight in estimating equations is closely related to the density power divergence (Basu et al., 1998) in which the density power divergence induces the density-weighted likelihood equations. However, the novelty of the proposed approach is that the density weight is considered within the framework of complete estimating equations rather than likelihood equations, which leads to tractable robust estimation algorithms for a variety of mixture models.

The remainder of this paper is organized as follows. In Section 2, we first describe the proposed WCE method for a general mixture model and derive a general EEE algorithm for solving the WCE. Then, the general algorithm is applied to the specific mixture models, a multivariate Gaussian mixture (Section 3), a mixture of experts (Section 4), and a multivariate skew normal mixture (Section 5). In Section 6, the proposed method is demonstrated for a multivariate Gaussian mixture and a skew normal mixture through the simulation studies. Section 7 illustrates practical advantages of the proposed method using real data. Finally, conclusions and discussions are presented in Section 8. The R code for implementing the proposed method is available at the GitHub repository (https://github.com/sshonosuke/RobMixture).

The Gaussian mixture models are the most famous and widely adopted mixture models for model-based clustering, even though the performance can be severely affected by outliers due to the light tails of the Gaussian (normal) distribution. Here, a new approach to the robust fitting of the Gaussian mixture model using the proposed weighted complete estimating equations is presented. Let $f({\text{\boldmath$y$}}_{i};{\text{\boldmath$\theta$}}_{k})=\phi_{p}({\text{\boldmath$y$}}_{i};{\text{\boldmath$\mu$}}_{k},{\text{\boldmath$\Sigma$}}_{k})$ denote the $p$-dimensional normal density $\phi_{p}$. It holds that

$$B({\text{\boldmath$\theta$}}_{k})=|{\text{\boldmath$\Sigma$}}_{k}|^{-\gamma/2}(2\pi)^{-p\gamma/2}(1+\gamma)^{-p/2},$$

and $C({\text{\boldmath$\theta$}}_{k})=(\textbf{0}_{p},{\rm vec}(C_{k2}))$, where $\textbf{0}_{p}$ is a $p$-dimensional vector of $0$ and

$\displaystyle C_{k2}=-\frac{\gamma}{2}|{\text{\boldmath$\Sigma$}}_{k}|^{-\gamma/2}{\text{\boldmath$\Sigma$}}_{k}^{-1}(2\pi)^{-p\gamma/2}(1+\gamma)^{-p/2-1}.$

Then, the weighted complete estimating equations for ${\text{\boldmath$\mu$}}_{k}$ and ${\text{\boldmath$\Sigma$}}_{k}$ are given by

	$\displaystyle\sum_{i=1}^{n}u_{ik}w({\text{\boldmath$y$}}_{i};{\text{\boldmath$\theta$}}_{k})({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k})=0$
	$\displaystyle\sum_{i=1}^{n}u_{ik}\Big{\{}w({\text{\boldmath$y$}}_{i};{\text{\boldmath$\theta$}}_{k}){\text{\boldmath$\Sigma$}}_{k}-w({\text{\boldmath$y$}}_{i};{\text{\boldmath$\theta$}}_{k})({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k})({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k})^{\top}\Big{\}}-g({\text{\boldmath$\Sigma$}}_{k})\Big{(}\sum_{i=1}^{n}u_{ik}\Big{)}{\text{\boldmath$\Sigma$}}_{k}=0,$

for $k=1,\ldots,K$, where $g({\text{\boldmath$\Sigma$}}_{k})=\gamma|{\text{\boldmath$\Sigma$}}_{k}|^{-\gamma/2}(2\pi)^{-p\gamma/2}(1+\gamma)^{-p/2-1}$. Hence, starting from some initial values ${\text{\boldmath$\Psi$}}^{(0)}=({\text{\boldmath$\pi$}}^{(0)},{\text{\boldmath$\theta$}}_{1}^{(0)},\ldots,{\text{\boldmath$\theta$}}_{K}^{(0)})$, the proposed EEE algorithm repeats the following two steps in the $s$th iteration:

• -

  E-step:   The standard E-step is left unchanged as

  $$u_{ik}^{(s)}=\frac{\pi_{k}^{(s)}\phi_{p}({\text{\boldmath$y$}}_{i};{\text{\boldmath$\mu$}}_{k}^{(s)},{\text{\boldmath$\Sigma$}}_{k}^{(s)})}{\sum_{\ell=1}^{K}\pi_{\ell}^{(s)}\phi_{p}({\text{\boldmath$y$}}_{i};{\text{\boldmath$\mu$}}_{\ell}^{(s)},{\text{\boldmath$\Sigma$}}_{\ell}^{(s)})}.$$

• -

  EE-step:   Compute the density weight $w_{ik}^{(s)}=w({\text{\boldmath$y$}}_{i};{\text{\boldmath$\theta$}}_{k}^{(s)})$ and then update the membership probabilities $\pi_{k}$’s as described in Section 2. The component-specific parameters $\theta_{k}$’s are updated as follows:

  	$\displaystyle{\text{\boldmath$\mu$}}_{k}^{(s+1)}=\frac{\sum_{i=1}^{n}u_{ik}^{(s)}w_{ik}^{(s)}{\text{\boldmath$y$}}_{i}}{\sum_{i=1}^{n}u_{ik}^{(s)}w_{ik}^{(s)}},$
  	$\displaystyle{\text{\boldmath$\Sigma$}}_{k}^{(s+1)}=\frac{\sum_{i=1}^{n}u_{ik}^{(s)}w_{ik}^{(s)}({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k}^{(s+1)})({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k}^{(s+1)})^{\top}}{\sum_{i=1}^{n}u_{ik}^{(s)}w_{ik}^{(s)}-g({\text{\boldmath$\Sigma$}}_{k}^{(s)})\sum_{i=1}^{n}u_{ik}^{(s)}},$

  and the mixing proportion $\pi_{k}$ is updated as (4).

Note that all the formulas in the above steps are obtained in closed form. This is one of the attractive features of the proposed method, as it does not require any computationally intensive methods such as Monte Carlo integration.

To choose reasonable starting values ${\text{\boldmath$\Psi$}}^{(0)}$, we employ existing robust estimation methods of Gaussian mixture models or clustering, such as trimmed robust clustering (TCL; Garcia-Escudero et al., 2008) and improper maximum likelihood (IML; Coretto and Hennig, 2016).

To stabilize the estimation of the variance-covariance matrix ${\text{\boldmath$\Sigma$}}_{k}$, it would be beneficial to impose the following eigen-ratio constraint:

$$\frac{\max_{j=1,\ldots,p}\max_{k=1,\ldots,K}\lambda_{j}({\text{\boldmath$\Sigma$}}_{k})}{\min_{j=1,\ldots,p}\min_{k=1,\ldots,K}\lambda_{j}({\text{\boldmath$\Sigma$}}_{k})}\leq c,$$

(6)

where $\lambda_{j}({\text{\boldmath$\Sigma$}}_{k})$ denotes the $j$th eigenvalue of the covariance matrix ${\text{\boldmath$\Sigma$}}_{k}$ in the $k$th component and $c$ is a fixed constant. When $c=1$, a spherical structure is imposed on ${\text{\boldmath$\Sigma$}}_{k}$, and a more flexible structure is allowed under a large value of $c$. To reflect the eigen-ratio constraint in the EEE algorithm, the eigenvalues of ${\text{\boldmath$\Sigma$}}_{k}^{(s)}$ can be simply replaced with the truncated version $\lambda_{j}^{\ast}({\text{\boldmath$\Sigma$}}_{k})$ where $\lambda_{j}^{\ast}({\text{\boldmath$\Sigma$}}_{k})=c$ if $\lambda_{j}({\text{\boldmath$\Sigma$}}_{k})>c$, $\lambda_{j}^{\ast}({\text{\boldmath$\Sigma$}}_{k})=\lambda_{j}({\text{\boldmath$\Sigma$}}_{k})$ if $c\theta_{c}\leq\lambda_{j}({\text{\boldmath$\Sigma$}}_{k})\leq c$ and $\lambda_{j}^{\ast}({\text{\boldmath$\Sigma$}}_{k})=c\theta_{c}$ if $\lambda_{j}({\text{\boldmath$\Sigma$}}_{k})<c{\text{\boldmath$\theta$}}_{c}$, where $\theta_{c}$ is an unknown constant that depends on $c$. The procedure of Fritz et al. (2013) is employed in this study. Note that if $c$ is set to a small value, the constraint for ${\text{\boldmath$\Sigma$}}_{k}$ is too severe such that the solution of the weighted estimating equation might not exist.

A mixture of experts model (Jacobs et al., 1991; McLachlan and Peel, 2000) is a useful tool for modelling nonlinear regression relationships. Models that include a simple mixture of normal regression models and their robust versions have been proposed in the literature (e.g. Bai et al., 2012; Song et al., 2014). In addition, the robust versions of a mixture of experts based on the non-normal distributions were considered by Nguyen and McLachlan (2016), and Chamroukhi (2016). The general model is described as

$$f_{M}(y_{i}|{\text{\boldmath$x$}}_{i};{\text{\boldmath$\Psi$}})=\sum_{k=1}^{K}g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{k})f(y_{i}|{\text{\boldmath$x$}}_{i};{\text{\boldmath$\theta$}}_{k}),$$

where ${\text{\boldmath$x$}}_{i}$ is a vector of covariates including an intercept term and $g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{k})$ is the mixing proportion such that $\sum_{k=1}^{K}g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{k})=1$. For the identifiability of the parameters, ${\text{\boldmath$\eta$}}_{K}$ is assumed to be an empty set, since $g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{K})$ is completely determined by ${\text{\boldmath$\eta$}}_{1},\ldots,{\text{\boldmath$\eta$}}_{K-1}$. A typical form of the continuous response variables adopts $f(y_{i}|{\text{\boldmath$x$}}_{i};{\text{\boldmath$\theta$}}_{k})=\phi(y_{i};{\text{\boldmath$x$}}_{i}^{\top}{\text{\boldmath$\beta$}}_{k},\sigma_{k}^{2})$ and $g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{k})=\exp({\text{\boldmath$x$}}_{i}^{\top}{\text{\boldmath$\eta$}}_{k})/\sum_{k=1}^{K}\exp({\text{\boldmath$x$}}_{i}^{\top}{\text{\boldmath$\eta$}}_{k})$. Let $\Psi$ be the set of the unknown parameters, ${\text{\boldmath$\eta$}}_{k}$, and ${\text{\boldmath$\theta$}}_{k}$. Compared to the standard mixture model (1), the mixing proportion $g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{k})$ is a function of ${\text{\boldmath$x$}}_{i}$ parameterized by ${\text{\boldmath$\eta$}}_{k}$.

As before, the latent variable $z_{i}$ such that $P(z_{i}=k)=g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{k})$ for $k=1,\ldots,K$ is introduced. Then, the complete weighted estimating equations are given by

$$\begin{split}&\sum_{i=1}^{n}u_{ik}\left\{w(y_{i};{\text{\boldmath$\theta$}}_{k})\frac{\partial}{\partial{\text{\boldmath$\theta$}}_{k}}\log f(y_{i};{\text{\boldmath$\theta$}}_{k})-C_{i}({\text{\boldmath$\theta$}}_{k})\right\}=0,\\ &\sum_{i=1}^{n}\frac{u_{ik}}{g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{k})}\frac{\partial g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{k})}{\partial{\text{\boldmath$\eta$}}_{k}}\frac{w(y_{i};{\text{\boldmath$\theta$}}_{k})}{B_{i}({\text{\boldmath$\theta$}}_{k})}-\sum_{i=1}^{n}\frac{u_{iK}}{g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{K})}\frac{\partial g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{k})}{\partial{\text{\boldmath$\eta$}}_{k}}\frac{w(y_{i};{\text{\boldmath$\theta$}}_{K})}{B_{i}({\text{\boldmath$\theta$}}_{K})}=0,\end{split}$$

where

	$\displaystyle C_{i}({\text{\boldmath$\theta$}}_{k})=\int f(t|{\text{\boldmath$x$}}_{i};{\text{\boldmath$\theta$}}_{k})w(t;{\text{\boldmath$\theta$}}_{k})\frac{\partial}{\partial{\text{\boldmath$\theta$}}_{k}}\log f(t|{\text{\boldmath$x$}}_{i};{\text{\boldmath$\theta$}}_{k})dt,$
	$\displaystyle\ \ \ \ B_{i}({\text{\boldmath$\theta$}}_{k})=\int f(t|{\text{\boldmath$x$}}_{i};{\text{\boldmath$\theta$}}_{k})w(t;{\text{\boldmath$\theta$}}_{k})dt.$

Starting from some initial values ${\text{\boldmath$\Psi$}}^{(0)}$, the proposed EEE algorithm repeats the following two steps in the $s$th iteration:

• -

  E-step:   The standard E-step is left unchanged as

  $$u_{ik}^{(s)}=\frac{g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{k}^{(s)})f(y_{i}|{\text{\boldmath$x$}}_{i};{\text{\boldmath$\theta$}}_{k}^{(s)})}{\sum_{\ell=1}^{K}g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{\ell}^{(s)})f(y_{i}|{\text{\boldmath$x$}}_{i};{\text{\boldmath$\theta$}}_{\ell}^{(s)})}.$$

• -

  EE-step:   Update ${\text{\boldmath$\eta$}}_{k}$ and ${\text{\boldmath$\theta$}}_{k}$ by solving the following equations:

  $$\begin{split}&\sum_{i=1}^{n}u_{ik}^{(s)}\left\{w(y_{i};{\text{\boldmath$\theta$}}_{k}^{(s)})\frac{\partial}{\partial{\text{\boldmath$\theta$}}_{k}}\log f(y_{i};{\text{\boldmath$\theta$}}_{k})-C_{i}({\text{\boldmath$\theta$}}_{k})\right\}=0,\\ &\sum_{i=1}^{n}\frac{u_{ik}^{(s)}}{g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{k})}\frac{\partial g}{\partial{\text{\boldmath$\eta$}}_{k}}\frac{w(y_{i};{\text{\boldmath$\theta$}}_{k}^{(s+1)})}{B_{i}({\text{\boldmath$\theta$}}_{k}^{(s+1)})}-\sum_{i=1}^{n}\frac{u_{iK}^{(s)}}{g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{K})}\frac{\partial g}{\partial{\text{\boldmath$\eta$}}_{k}}\frac{w(y_{i};{\text{\boldmath$\theta$}}_{K}^{(s+1)})}{B_{i}({\text{\boldmath$\theta$}}_{K}^{(s+1)})}=0.\end{split}$$

  (7)

When the mixture components are the normal linear regression models given by $f(y_{i}|{\text{\boldmath$x$}}_{i};{\text{\boldmath$\theta$}}_{k})=\phi(y_{i};{\text{\boldmath$x$}}_{i}^{\top}{\text{\boldmath$\beta$}}_{k},\sigma_{k}^{2})$ with ${\text{\boldmath$\theta$}}_{k}=({\text{\boldmath$\beta$}}_{k}^{\top},\sigma_{k}^{2})$, the bias correction terms $C_{i}({\text{\boldmath$\theta$}}_{k})\equiv(\textbf{0}_{p},C_{i2}({\text{\boldmath$\theta$}}_{k}))$ and $B_{i}({\text{\boldmath$\theta$}}_{k})$ can be analytically obtained as

$$C_{i2}({\text{\boldmath$\theta$}}_{k})=-(\gamma/2)(\sigma_{k}^{2})^{-1-\gamma/2}(2\pi)^{-\gamma/2}(1+\gamma)^{-3/2}$$

$$B_{i}({\text{\boldmath$\theta$}}_{k})=(2\pi\sigma_{k}^{2})^{-\gamma/2}(1+\gamma)^{-1/2}.$$

The first equation in (7) can be solved analytically, and the closed-form updating steps similar to those in Section 3 can be obtained. On the other hand, the second equation is a function of ${\text{\boldmath$\eta$}}_{1},\ldots,{\text{\boldmath$\eta$}}_{K-1}$ and cannot be solved analytically. Note that the solution corresponds to the maximizer of the weighted log-likelihood function of the multinomial distribution given by

$$\sum_{i=1}^{n}\sum_{k=1}^{K}u_{ik}^{(s)}\frac{w({\text{\boldmath$y$}}_{i};{\text{\boldmath$\theta$}}_{k}^{(s+1)})}{B_{i}({\text{\boldmath$\theta$}}_{k}^{(s+1)})}\log g({\text{\boldmath$x$}}_{i};{\text{\boldmath$\eta$}}_{k}),$$

since its first-order partial derivatives with respect to ${\text{\boldmath$\eta$}}_{k}$ are reduced to the second equation in (7). Thus, the updating step for ${\text{\boldmath$\eta$}}_{k}$ can be readily carried out.

In setting the initial values of the EEE algorithm, we first randomly split the data into $K$ groups and apply some existing robust methods to estimate ${\text{\boldmath$\theta$}}_{k}$ for $k=1,\ldots,K$, which can be adopted for the initial values ${\text{\boldmath$\theta$}}_{k}^{(0)}$ of ${\text{\boldmath$\theta$}}_{k}$. For example, when $f(\cdot;{\text{\boldmath$\theta$}}_{k})$ is the normal linear regression as adopted in Section 7.1, we employ an M-estimator available in the rlm function in R package ‘MASS’ (Venables and Ripley, 2002). For the initial values of ${\text{\boldmath$\eta$}}_{k}$, we suggest simply using ${\text{\boldmath$\eta$}}_{k}^{(0)}=(0,\ldots,0)$.

We next consider the use of the $p$-dimensional skew normal distribution (Azzalini and Valle, 1996) for the $k$th component. The mixture model based on skew normal distributions is more flexible than a multivariate Gaussian mixture, especially when the cluster-specific distributions are not symmetric but skewed. Several works have been conducted on the maximum likelihood estimation of a skew normal mixture (Lin et al., 2007; Lin, 2009). Despite the flexibility in terms of skewness, because the skew normal distribution still has light tails as the normal distribution, the skew normal mixture is also sensitive to outliers. Alternative mixture models using heavy-tailed and skewed distributions have also been proposed in the literature (e.g. Lee and McLachlan, 2017; Morris et al., 2019). Here, we consider the robust fitting of the skew normal mixture within the proposed framework.

We first note that the direct application of the skew normal distribution to the proposed EEE algorithm would be computationally intensive since the bias correction term $C({\text{\boldmath$\theta$}}_{k})$ cannot be obtained analytically, and a Monte Carlo approximation would be required in each iteration. Instead, we provide a novel algorithm that exploits the stochastic representation of the multivariate skew normal distribution. The resulting algorithm is a slight extension of the proposed algorithm in Section 2. Früwirth-Schnatter and Pyne (2010) provided the following hierarchical representation for ${\text{\boldmath$y$}}_{i}$ given $z_{i}$:

$$\begin{split}&{\text{\boldmath$y$}}_{i}|(z_{i}=k)={\text{\boldmath$\mu$}}_{k}+{\text{\boldmath$\psi$}}_{k}v_{ik}+{\varepsilon}_{ik},\\ &\epsilon_{ik}\sim N(0,{\text{\boldmath$\Sigma$}}_{k}),\quad v_{ik}\sim N^{+}(0,1),\end{split}$$

(8)

where ${\text{\boldmath$\mu$}}_{k}$ is the $p\times 1$ vector of location parameters, ${\text{\boldmath$\psi$}}_{k}$ is the $p\times 1$ vector of the skewness parameters, and ${\text{\boldmath$\Sigma$}}_{k}$ is the covariance matrix. Here, $N^{+}(a,b)$ denotes the truncated normal distribution on the positive real line with mean and variance parameters $a$ and $b$, respectively, and its the density function $\phi^{+}(x;a,b^{2})$ is given by

$\displaystyle\phi^{+}(x;a,b)=\frac{1}{\Phi(a/b)}\phi(x;a,b)I(x\geq 0),$

where $\Phi(\cdot)$ is the distribution function of the standard normal distribution. By defining

$${\text{\boldmath$\Omega$}}_{k}={\text{\boldmath$\Sigma$}}_{k}+{\text{\boldmath$\psi$}}_{k}{\text{\boldmath$\psi$}}_{k}^{\top},\quad{\text{\boldmath$\alpha$}}_{k}=\frac{{\text{\boldmath$\omega$}}_{k}{\text{\boldmath$\Omega$}}_{k}^{-1}{\text{\boldmath$\psi$}}_{k}}{\sqrt{1-{\text{\boldmath$\psi$}}_{k}^{\top}{\text{\boldmath$\Omega$}}_{k}^{-1}{\text{\boldmath$\psi$}}_{k}}},$$

where ${\text{\boldmath$\omega$}}_{k}=({\text{\boldmath$\Omega$}}_{11}^{1/2},\dots,{\text{\boldmath$\Omega$}}_{pp}^{1/2})$. The probability density function of the multivariate skew normal distribution of Azzalini and Valle (1996) is given by

$$f_{SN}({\text{\boldmath$y$}};{\text{\boldmath$\mu$}}_{k},{\text{\boldmath$\Omega$}}_{k},{\text{\boldmath$\alpha$}}_{k})=2\phi_{p}({\text{\boldmath$y$}};{\text{\boldmath$\mu$}}_{k},{\text{\boldmath$\Omega$}}_{k})\Phi({\text{\boldmath$\alpha$}}_{k}^{\top}{\text{\boldmath$\omega$}}^{-1}_{k}({\text{\boldmath$y$}}-{\text{\boldmath$\mu$}}_{k})).$$

(9)

From (8), the conditional distribution of ${\text{\boldmath$y$}}_{i}$ given both $z_{i}$ and $v_{ik}$ is normal, namely, ${\text{\boldmath$y$}}_{i}|(z_{i}=k),v_{ik}\sim N({\text{\boldmath$\mu$}}_{k}+{\text{\boldmath$\psi$}}_{k}v_{ik},{\text{\boldmath$\Sigma$}}_{k})$, so that the bias correction terms under given $z_{i}$ and $v_{ik}$ can be easily obtained in the same way as in the case of the Gaussian mixture models. Therefore, we consider the following complete estimating equations for the parameters conditional on both $z_{i}$ and $v_{ik}$:

$$\begin{split}&\sum_{i=1}^{n}u_{ik}w_{ik}({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k}-{\text{\boldmath$\psi$}}_{k}v_{ik})=0,\\ &\sum_{i=1}^{n}u_{ik}w_{ik}\Big{\{}{\text{\boldmath$\Sigma$}}_{k}-({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k}-{\text{\boldmath$\psi$}}_{k}v_{ik})({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k}-{\text{\boldmath$\psi$}}_{k}v_{ik})^{\top}\Big{\}}-g({\text{\boldmath$\Sigma$}}_{k})\left(\sum_{i=1}^{n}u_{ik}\right){\text{\boldmath$\Sigma$}}_{k}=0,\\ &\sum_{i=1}^{n}u_{ik}v_{ik}w_{ik}({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k}-{\text{\boldmath$\psi$}}_{k}v_{ik})=0,\\ \end{split}$$

where $w_{ik}\equiv w({\text{\boldmath$y$}}_{i};{\text{\boldmath$\theta$}}_{k})$ and $g({\text{\boldmath$\Sigma$}}_{k})$ have the same forms as in the Gaussian mixture case.

The proposed EEE algorithm for the robust fitting of the skew normal mixture is obtained after some modification to the normal case. Since the complete estimating equations contain the additional latent variables $v_{ik}$, some additional steps are included to impute the moments of $v_{ik}$ in the equations. The moments are those of the conditional posterior distribution of $v_{ik}$, given $z_{i}=k$, which is $N^{+}(\delta_{ik},\tau_{ik}^{2})$, where

$$\delta_{ik}=\frac{{\text{\boldmath$\psi$}}_{k}^{\top}{\text{\boldmath$\Sigma$}}_{k}^{-1}({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k})}{{\text{\boldmath$\psi$}}_{k}^{\top}{\text{\boldmath$\Sigma$}}_{k}^{-1}{\text{\boldmath$\psi$}}_{k}+1},\ \ \ \ \tau_{ik}^{2}=\frac{1}{{\text{\boldmath$\psi$}}_{k}^{\top}{\text{\boldmath$\Sigma$}}_{k}^{-1}{\text{\boldmath$\psi$}}_{k}+1}.$$

We define the following quantities:

$$\begin{split}&t^{2}_{ik}=\left(\gamma{\text{\boldmath$\psi$}}_{k}^{\top}{\text{\boldmath$\Sigma$}}_{k}^{-1}{\text{\boldmath$\psi$}}_{k}+\frac{1}{\tau^{2}_{ik}}\right)^{-1},\ \ \ \ m_{ik}=t^{2}_{ik}\left\{\gamma{\text{\boldmath$\psi$}}_{k}^{\top}{\text{\boldmath$\Sigma$}}_{k}^{-1}({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k})+\frac{\delta_{ik}}{\tau_{ik}^{2}}\right\},\\ &U_{ik}=|{\text{\boldmath$\Sigma$}}_{k}|^{-\gamma/2}(2\pi)^{-\gamma p/2}\frac{\Phi(m_{ik}/t_{ik})}{\Phi(\delta_{ik}/\tau_{ik})}\left(\frac{t_{ik}^{2}}{\tau_{ik}^{2}}\right)^{1/2}\\ \ &\times\exp\left\{-\frac{\gamma}{2}({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k})^{\top}{\text{\boldmath$\Sigma$}}_{k}^{-1}({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k})-\frac{\delta^{2}_{ik}}{2\tau^{2}_{ik}}+\frac{m_{ik}^{2}}{2t_{ik}^{2}}\right\}.\end{split}$$

(10)

Note that we have for $j=0,1$ and $2$

$\displaystyle E[u_{ik}w_{ik}v_{ik}^{j}|{\text{\boldmath$y$}}_{i}]$

$\displaystyle=E[E[u_{ik}w_{ik}v_{ik}^{j}|{\text{\boldmath$y$}}_{i},u_{ik}]|{\text{\boldmath$y$}}_{i}]=E[u_{ik}|{\text{\boldmath$y$}}_{i}]E[w_{ik}v_{ik}^{j}|u_{ik},{\text{\boldmath$y$}}_{i}].$

Since the conditional distribution of $v_{ik}$ given $(u_{ik},{\text{\boldmath$y$}}_{i})$ is $N^{+}(\delta_{ik},\tau_{ik}^{2})$, we have

	$\displaystyle E[v_{ik}^{j}w_{ik}|u_{ik},{\text{\boldmath$y$}}_{i}]$
	$\displaystyle=\int_{0}^{\infty}|{\text{\boldmath$\Sigma$}}_{k}|^{-\gamma/2}(2\pi)^{-\gamma p/2}\exp\left\{-\frac{\gamma}{2}({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k}-v_{ik}{\text{\boldmath$\psi$}}_{k})^{\top}{\text{\boldmath$\Sigma$}}_{k}^{-1}({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k}-v_{ik}{\text{\boldmath$\psi$}}_{k})\right\}$
	$\displaystyle\ \ \ \ \ \times\frac{v_{ik}^{j}}{\Phi(\delta_{ik}/\tau_{ik})}(2\pi\tau_{ik}^{2})^{-1/2}\exp\left\{-\frac{(v_{ik}-\delta_{ik})^{2}}{2\tau_{ik}^{2}}\right\}dv_{ik}$
	$\displaystyle=|{\text{\boldmath$\Sigma$}}_{k}|^{-\gamma/2}(2\pi)^{-\gamma p/2}\frac{\Phi(m_{ik}/t_{ik})}{\Phi(\delta_{ik}/\tau_{ik})}\left(\frac{t_{ik}^{2}}{\tau_{ik}^{2}}\right)^{1/2}\int_{0}^{\infty}v_{ik}^{j}\phi_{+}(v_{ik};m_{ik},t_{ik}^{2})dv_{ik}$
	$\displaystyle\ \ \ \ \ \ \times\exp\left\{-\frac{\gamma}{2}({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k})^{\top}{\text{\boldmath$\Sigma$}}_{k}^{-1}({\text{\boldmath$y$}}_{i}-{\text{\boldmath$\mu$}}_{k})-\frac{\delta^{2}_{ik}}{2\tau^{2}_{ik}}+\frac{m_{ik}^{2}}{2t_{ik}^{2}}\right\},$

and it follows from Lin et al. (2007) that

	$\displaystyle\int_{0}^{\infty}v_{ik}^{1}\phi_{+}(v_{ik};m_{ik},t_{ik}^{2})dv_{ik}=m_{ik}+t_{ik}\frac{\phi(m_{ik}/t_{ik})}{\Phi(m_{ik}/t_{ik})},$
	$\displaystyle\int_{0}^{\infty}v_{ik}^{2}\phi_{+}(v_{ik};m_{ik},t_{ik}^{2})dv_{ik}=m_{ik}^{2}+t_{ik}^{2}+m_{ik}t_{ik}\frac{\phi(m_{ik}/t_{ik})}{\Phi(m_{ik}/t_{ik})}.$

These expressions are applied to the analytical updating steps in E-step. Starting from some initial values of the parameters, ${\text{\boldmath$\Psi$}}^{(0)}=({\text{\boldmath$\pi$}}^{(0)},{\text{\boldmath$\theta$}}_{1}^{(0)},\ldots,{\text{\boldmath$\theta$}}_{K}^{(0)})$, the proposed EEE algorithm iteratively updates the parameters in the $s$th iteration as follows: