{'fraction_non_alphanumeric': 0.04412397144413544, 'fraction_numerical': 0.013131110410275523, 'mean_word_length': 4.496375669713205, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Real Elliptically Skewed Distributions and Their Application to Robust Cluster Analysis

This article proposes a new class of Real Elliptically Skewed (RESK) distributions and associated clustering algorithms that allow for integrating robustness and skewness into a single unified cluster analysis framework. Non-symmetrically distributed and heavy-tailed data clusters have been reported in a variety of real-world applications. Robustness is essential because a few outlying observations can severely obscure the cluster structure. The RESK distributions are a generalization of the Real Elliptically Symmetric (RES) distributions. To estimate the cluster parameters and memberships, we derive an expectation maximization (EM) algorithm for arbitrary RESK distributions. Special attention is given to a new robust skew-Huber M-estimator, which is also the maximum likelihood estimator (MLE) for the skew-Huber distribution that belongs to the RESK class. Numerical experiments on simulated and real-world data confirm the usefulness of the proposed methods for skewed and heavy-tailed data sets.

Abstract-This article proposes a new class of Real Elliptically Skewed (RESK) distributions and associated clustering algorithms that integrate robustness and skewness into a single unified cluster analysis framework. Non-symmetrically distributed and heavy-tailed data clusters have been reported in a variety of real-world applications. Robustness is essential because a few outlying observations can severely obscure the cluster structure. The RESK distributions are a generalization of the Real Elliptically Symmetric (RES) distributions. To estimate the cluster parameters and memberships, we derive an expectation maximization (EM) algorithm for arbitrary RESK distributions. Special attention is given to a new robust skew-Huber M-estimator, which is also the approximate maximum likelihood estimator (MLE) for the skew-Huber distribution, that belongs to the RESK class. Numerical experiments on simulated and real-world data confirm the usefulness of the proposed methods for skewed and heavy-tailed data sets.

I. INTRODUCTION
F INITE mixture models have been extensively used in cluster analysis to estimate the probability density function (pdf) of a given data set [1]- [5]. It is well-known that estimators that have been derived under a Gaussian data assumption are strongly affected by outliers, which frequently occur in real-world applications [6]- [8]. A popular remedy is to choose a heavy-tailed finite mixture model, such as, a member of the real elliptically symmetric (RES) family of distributions. The associated maximum-likelihood-estimators (MLE) can handle outliers and heavy-tailed noise. M-estimators provide further robustness, as they decouple the estimation from a specific distribution. The underlying idea is to replace the negative loglikelihood function by a robustness-inducing objective function [9]- [12].
While the RES distributions offer a high degree of flexibility in modeling different types of heavy-tailedness in the data distribution, they still depend on the assumption that the data is symmetrically distributed. However, in medical [13], [14], localization [15], [16] or financial [17] applications, it is often the case that the data is heavy-tailed and nonsymmetrically distributed. The focus of this work, therefore, lies on mixture models that account for such cases. The proposed robust clustering algorithms are designed to provide reliable results, even if the data distribution is skewed, heavytailed, and if the data set contains outliers. A widely used C. A. Schroth  estimation algorithm for mixture models is the expectation maximization (EM) algorithm [1]. The standard EM algorithm is based on a Gaussian assumption [18], hence, it is not robust against outliers. Some efforts were made to robustify the EM by outlier detection and removal [19]- [22], assuming t-distributed/heavy-tailed data [23]- [25], using an additional component in a mixture model [26], [27], including contaminated Gaussian distributions [28] or robustifying the estimation of the cluster centers and covariances [29]. Recently, an EM algorithm has been derived that is applicable to any RES distribution [30]. An important member of the RES family, is the so-called Huber distribution, for which the popular Huber M-estimator is the MLE. Multivariate skewed distributions include the skew-Gaussian distribution [31], [32], the skew-t distribution [33] and skew-Laplace distribution [34]. A summary of different skewed distributions can be found in [35], [36]. In contrast to RES distributions, there are multiple formulations for skewed distributions and, generally speaking, none of the specific formulations seems to be preferable over the others [37], [38]. In this work, we use the formulation introduced by [39] because it allows for deriving an EM algorithm with explicit solutions for the model parameters. This does not seem to be possible, for example, for the generalized multivariate skewed distributions introduced in [40], [41].
EM algorithms are available for a wide range of specific skewed distributions, e.g., for the skew-Gaussian [39], [42]- [47], the skew-t [45], [46], [48], [49], the skew power exponential [50] and the generalized hyperbolic [51]. Nevertheless to the best of our knowledge, a unified estimation framework for a broad class of skewed distributions has not yet been developed.
To fit a finite mixture model to a given data set, the EM algorithm requires the number of clusters to be known, or estimated. A popular strategy is to use model selection criteria, such as, the Bayesian Information Criterion (BIC) derived by Schwarz [52], [53]. Alternatively, a Bayesian cluster enumeration approach that selects the model that is maximum a posteriori most probable [54], [55] may be beneficial, because it accounts for the structure of the clustering problem in the penalty term.
Based on [39], [45]- [47], our main contribution is to propose a formulation for Real Elliptically Skewed (RESK) distributions and to derive an EM algorithm for such distributions. The RESK family of distributions is a generalization of the RES family of distributions, which are included within the RESK class by choosing an appropriate value for the skewness factor. For robustness purposes, our framework is then generalized to include asymmetric M-estimators, where arXiv:2006.16671v2 [eess.SP] 21 Jun 2021 special attention is given to a new skew-Huber's M-estimator. Numerical experiments demonstrate increased parameter estimation and cluster enumeration performance for skewed and heavy-tailed data.
The paper is organized as follows. Section II briefly revisits RES distributions and loss functions. Section III introduces the proposed RESK family of distributions with special attention given to the skew-Huber distribution. Section IV derives the proposed EM algorithm for RESK distributions and Section V briefly describes cluster enumeration for RESK distributions based on the Schwarz BIC. Simulations and real-world examples are provided and discussed in Section VI. Finally, conclusions are drawn and some possible future steps are mentioned in Section VII.
Notation: Normal-font letters (a, A) denote a scalar, bold lowercase (a) a vector and bold uppercase (A) a matrix; calligraphic letters (A) denote a set, with the exception of L, which is used for the likelihood function; R is the set of real numbers and R r×1 , R r×r the set of column vectors of size r×1 and matrices of size r×r, respectively; A −1 is the matrix inverse; A is the matrix transpose; |a| is the absolute value of a scalar; |A| is the determinant of a matrix; ⊗ represents the Kronecker product; vec(·) is the vectorization operator, D is the duplication matrix and vech(·) is the vector half operator as defined in [56]. The vech operator takes a symmetric r × r matrix and stacks the lower/upper triangular half into a single vector of length r 2 (r + 1).

A. RES Distributions
Assuming that the observed data x ∈ R r×1 follows a RES distribution with centroid µ ∈ R r×1 and positive definite symmetric scatter matrix S ∈ R r×r , the pdf of x is given by where t = (x − µ) S −1 (x − µ) is the squared Mahalanobis distance, see [7, p. 109] and [41]. The function g, often referred to as the density generator, is defined by where h(t; r) is a function for which holds. Note that h(t; r) can be a function of multiple parameters, not only of r.

B. Loss Functions
Assuming an observation of N independent and identically distributed (iid) random variables denoted by X = {x 1 , . . . , x N }, the associated likelihood function is given by with t n = (x n − µ) S −1 (x n − µ). The MLE minimizes the negative log-likelihood function where the ML loss function [7, p. 109] is defined as The corresponding first and second derivatives are denoted, respectively, by The key idea of M-estimation [9] is to replace the ML loss function ρ ML (t n ) in Eq. (6) with a more general loss function ρ(t n ). The aim of robust M-estimators is not necessarily to be optimal for a single specific distribution, but to provide near-optimal performance within a neighborhood of distributions including heavy-tailed distributions. The weight function ψ(t n ) can be designed according to desired robustness properties. For example, the very popular Huber weight function down weights outlying data points while giving full weight to points with a small Mahalanobis distance. Interestingly, Huber's estimator is not only a popular M-estimator, but it also has been shown to be the MLE for a specific RES distribution. An overview of some exemplary loss functions and their derivatives can be found in Tables I and II. Since the Gaussian and t distribution are well-known they are not further discussed, but for the Huber distribution a brief discussion is provided in Section III-C.

III. RESK DISTRIBUTIONS A. Definition of RESK Distributions
This section proposes the class of RESK distributions as a generalization of RES distributions to allow for skewness. Well-known skewed distributions, such as, the skew-Gaussian and the skew-t [14], [41], [45]- [47] distribution are included as special cases. The aim is to provide a unified framework to derive robust clustering algorithms for skewed data. Based on the RESK distributions, the user is given a large flexibility in designing algorithms that fit to skewed and heavy-tailed data. This is because the proposed formulation allows for the derivation of an explicit solution for the EM algorithm for arbitrary RESK distributions. As an example, a new clustering algorithm by skewing Huber's distribution is derived.
Assuming that the observed data x ∈ R r×1 follows a RESK distribution, let ξ ∈ R r×1 be the 'centroid', λ ∈ R r×1 the skewness factor and let S ∈ R r×r be the positive definite symmetric scatter matrix. Then, it is possible to define the RESK distribution as a RES distribution multiplied with its univariate cdf with the skewed scatter matrix the skewed squared Mahalanobis distance and with ψ (t) from Table I. The normalization factor '2' in Eq. (8) arises naturally from the multiplication with the cdf of a RES distribution, because every cdf of a RES distribution has a mean value of 1/2. This result can be easily verified by taking the symmetry of the RES distributions into account. Hence, the same result is obtained by evaluating F (0) = 1/2. To simplify the notation, the argument of κ (t) is dropped and the following scalar variables are introduced: Figure 1a illustrates the flexibility in modeling skewed data by plotting the univariate skew-Gaussian distribution for different values of the skewness parameter λ. In the multivariate case, the skewness vector λ allows for different amounts of skewness in every dimension. Figure 1b shows the influence of the centroid ξ, which for the univariate case simply shifts the distribution to the left or right. In the multivariate case, ξ shifts independently the distribution in each dimension. While the figure illustrates these parameters for the skew-Gaussian distribution, their meaning is identical for other RESK distributions. We next briefly revisit some known special cases of RESK distributions before introducing the skew-Huber distribution. All required functions are summarized in Tables I and II. B. Examples of Known Existing RESK Distributions 1) RES Distributions: The RESK distribution has the property that, for λ = 0, it is equal to a RES distribution as defined in Eq. (1), because η = 0 , τ = 1 and κ (t) = 0, which leads to F (0) = 1/2. Also, one can find that ξ = µ, hence the RESK distribution is a generalization of the RES distribution. This can be also observed in Figure 1a. By comparing the black dashed curve with the red curve, it becomes clear that, for λ = 0, the skew-Gaussian is equal to the Gaussian distribution.
2) Skew-Gaussian Distribution: The pdf is defined by with φ(x|µ, Σ) being the multivariate pdf and Φ(x) being the univariate cdf of the Gaussian distribution.
3) Skew-t Distribution: The pdf is defined by where t(x|µ, Ψ, ν) is the multivariate pdf and T ν (x) is the univariate cdf of a t distribution with associated degree of freedom parameter ν. During our work with the skew-t distribution we came across some inconsistencies regarding the definition of t and η τ of the skew-t distribution between [46] and [47]. These inconsistencies are corrected here and the correct expressions are given in Eqs. (10) and (12).

C. Proposed Skew-Huber Distribution
Based on the definition in Eq. (8), the skew-Huber distribution can now be introduced as with some tuning constant c ≥ 0 and To obtain Fisher consistency where F χ 2 r (·) is the Chi-square cumulative distribution function with degree of freedom r. For h(t; r, c) to be a valid pdf the normalization factor, according to [41], has to be calculated as with the gamma function Γ(·), the upper incomplete gamma function Γ(·, ·) and Euler's number e. The density generator, see Eq.
(1), of a Huber distribution now follows as with (21) Then, the multivariate Huber pdf is obtained as 2) Univariate cdf of Huber Distribution: We start with the univariate pdf and integrate it. This leads to with the constraint c 2 2b > 1 and When deriving estimators based on the Huber distribution, [7, p. 116] suggest to choose c 2 as the q H th upper quantile of a χ 2 r distribution The value of q H is chosen to trade-off efficiency under a Gaussian model and robustness against outliers. We can proceed similarly for the skew-Huber, but from c 2 2b > 1 in Eq. (25), we have the constraint that q H > 0.703. In Figure 1c, the skew-Huber distribution is plotted over different values of q H from Eq. (27). For lower values of q H , probability mass is shifted to the tails of the distribution. As q H → 1 and c → ∞ the skew-Huber becomes equal to the skew-Gaussian distribution. Additionally, in Figure 1d the influence of the weight function ψ(t) is shown. Similarly, as q H → 1 and c → ∞ the weight function becomes equal to the Gaussian weight function.

IV. PROPOSED EXPECTATION MAXIMIZATION ALGORITHM FOR A MIXTURE OF RESK DISTRIBUTIONS
In this section, we derive an EM algorithm to estimate the RESK mixture model parameters and the cluster memberships of the data vectors x n by alternating the E-and M-steps until convergence. The proposed approach is summarized in Algorithm 1, where we provide a unified EM clustering framework that is beneficial for skewed and heavy-tailed data sets. The user simply chooses a RESK density generator g(t) and determines the associated weight function ψ(t). We have implemented some specific variants of the EM algorithm, for example the skew-Huber, based on the values from Tables I  and II. However, the framework is not limited to these cases, and therefore it is very flexible in modeling different data distributions.
The parameter vector, for each cluster m = 1, . . . , l, is * The calculation of Ψ(x) can be numerically unstable, due to the fact that Φ(x) tends to zero for x < −37.
An approximation for these values can be found in [46, p. 29].
l is the number of clusters and q = r 2 (r + 5) is the number of parameters per cluster. Because S m is symmetric, it has only r 2 (r+1) unique elements. Therefore, vech(S m ) has to be used in the estimation step [58, p. 367]. Accordingly, we derive an EM algorithm for RESK distributions in the following.
Based on Equation (8), the mixture of l RESK distributions is formulated as and the incomplete data log-likelihood follows as with γ m being the mixing coefficient and Following the derivation of the EM algorithm from [18], a latent binary variable z = [z 1 , . . . , z l ] is introduced, which has a single component equal to 1 and all other components equal to 0. The conditional distribution of x given z, becomes where the prior of the latent variables is given by and the joint distribution equals Now, the complete data log-likelihood, with X and Z = {z 1 , . . . , z N }, is formulated as where z nm denotes the mth element of z n . For the E-step the expected value of the complete data log-likelihood with respect to the latent variables is calculated as The calculation of the ML estimates when using the EM algorithm requires fixed values for the variable v nm and additionally for e 0,nm , e 1,nm and e 2,nm from Eqs. (53), (54) and (64). In the E-step of the EM algorithm v nm is fixed with respect to θ m by calculating the expected value of z nm . A latent variable approach would be also favorable for e 0,nm , e 1,nm and e 2,nm , but it turns out, that for general RESK distributions this is not feasible. Hence, e 0,nm , e 1,nm and e 2,nm are approximated by the estimated valuesθ m of the previous iteration of the EM algorithm into Eqs. (53), (54) and (64) which yields the approximationsẽ 0,nm ,ẽ 1,nm andẽ 2,nm . If the approximations hold exactly, which is the case e.g. for the skew-Gaussian distribution, the ML estimate is obtained.
Otherwise an approximate ML estimate is obtained. In the context of M-estimation, where the primary aim is trading off near-optimality at a nominal distribution for robustness over a family of distributions [7], [9], such an approximation is well motivated and of high practical value. Based on the approximations, the expected value of the complete data log-likelihood, is maximized, which results in the following estimates for the M-step: A detailed derivation of the M-step is given in Appendix A.
V. CLUSTER ENUMERATION This section gives a brief overview of the Schwarz BIC [52], [53], which we use to estimate the number of clusters. There exist alternative criteria, e.g. [30], [59], [60], but in this paper the focus does not lie on cluster enumeration. Hence, only the log-likelihood function is adapted to RESK distributions and the Schwarz BIC is applied, i.e.
with M l being a candidate model with l ∈ {L min , . . . , L max } clusters. The true number of clusters K is assumed to lie within L min ≤ K ≤ L max . For the RES case, the loglikelihood function is Algorithm 1: EM algorithm for RESK distribution Input: X , i max , l, g(t), ψ(t) Output:ξ m ,Ŝ m ,λ m ,γ m for m = 1, . . . , l do Initializeξ (0) m , e.g. with K-means or K-medoids with Eqs. (35), (53), (54) and (64) Calculate log-likelihood: For the RESK case, the log-likelihood function is derived as The number of clusters is estimated by evaluatinĝ

VI. EXPERIMENTAL RESULTS
This section reports on our numerical experiments to evaluate the proposed EM algorithms based on RESK distributions in terms of (i) their parameter estimation and clustering performance for skewed and heavy-tailed mixture models; (ii) their robustness against outliers; (iii) their convergence speed in terms of the number of required iterations; and (iv) their ability to model real-world data. We compare the newly proposed skew-Huber estimator to other known estimators, such as the skew-Gaussian, skew-t [46], and the classic Gaussian, t-and Huber estimators. An implementation of the proposed EM algorithm and the conducted simulations is available at: https: //github.com/schrchr/Real-Elliptically-Skewed-Distributions.

A. Dataset 1
The first data set is comprised of K = 3 skew-Gaussian distributed clusters. It is generated with the algorithm from [46, p. 23], and x k ∼ N s (ξ k , Σ k , λ k ), k = 1, 2, 3. The number of data points per cluster is specified as An exemplary realization is given in Figure 2a.

C. Outlier Generation
To study robustness against outliers, an -contamination model is used, where denotes the percentage of randomly replaced data points. To make the clustering problem welldefined, it is assumed that the outliers do not form a cluster. Therefore, a uniform distribution is chosen as the outlier generating distribution because it produces isolated data points that are placed with equal probability on a sufficiently large area compared to the clustered data. For Dataset 1 and Dataset 2, the outliers are uniformly distributed in the range of [-15, 45], [-20, 30] in x, y-dimension.

D. Choice of Tuning Parameters
Based on [7, p. 121], for Huber's and skew-Huber's estimator, a value of q H = 0.8 is chosen for all simulations, which leads to c = 1.282. For all t and skew-t distributions, ν = 3 is selected for the degree of freedom. Both are common choices in M-estimation to trade-off robustness against outliers and efficiency for Gaussian data. To evaluate the choice of ν = 3, we also report on an oracle version, where the degree of freedom is selected from the set ν = {1, 2, 3, 4, 5, 7, 9} so that it provides the optimal performance. E. Performance Measures 1) Kullback-Leibler Divergence: To evaluate the parameter estimation performance of the proposed method, we use the Kullback-Leibler (KL) divergence, which measures the difference between two pdfs [61, p. 34]. It is defined as with p(x i ) representing the pdf that generates the outlierfree data and q(x i ) representing the pdf estimated by the EM algorithm. Even if there are outliers in the data, the aim of robust methods is to provide information about p(x i ). Therefore, the KL divergence between p(x i ) and q(x i ) is computed.
2) Confusion Matrix: To evaluate the correctness in labeling the non-outlying data points, the confusion matrix for these points is reported. 3) Sensitivity Curve: The sensitivity curve (also known as empirical influence function [7]) shows the effect of a single outlier as a function of its position in the feature space. The concept of looking at single outliers is related to infinitesimal robustness, or qualitative robustness, as measured by the influence function. A robust method is expected to provide a nearly optimal performance and is not largely influenced by a single outlier, irrespective of the value that the outlier takes. In our experiments, the sensitivity curves are calculated as follows. One randomly chosen sample is replaced with an outlier from the range [-15, 45], [-20, 30] in x, y-dimension. For each point of a grid of outlier values, the proposed EM algorithm is executed. The resulting estimated pdf is then compared with the true (outlier-free) generating pdf by calculating the KL divergence, which is averaged over 200 Monte Carlo iterations with N K = 50. 4) Breakdown Behavior: Showing robustness against a single outlier is not sufficient in cluster analysis, where we may have positive a fraction of data points that do not follow the cluster structure of the remaining data points. This is several experiments with varying amounts of outliers are conducted. To analyze robustness against a fraction of outliers, data is generated from a mixture model with a nominal distribution (e.g. the Gaussian or the t) and a contaminating (outlier generating) distribution. The level of contamination is increased until all methods break down. iterations. In the first row, the number of clusters is estimated with the Schwarz BIC from Eq. (40) and in the second row the true number of clusters is assumed to be known. By comparing the first column (RES model) with the second column (RESK model), it becomes clear that, for this simulation, the proposed skew-Huber model performs better than all other RES and RESK models. This is the case, because of Huber's weight function ψ(t n ) gives full weight to the data points with a squared distance smaller than c 2 while heavily downweighting the outliers. One can also note that the skew-t with ν = 3 is outperformed in all cases, and even the oracle method skew-t, opt is outperformed for almost all cases when K is assumed to be known. Until approximately 3% of outliers, the KL divergences for the true and estimated number of clusters are almost constant and similar, because in this range the Schwarz BIC has a very high empirical probability to detect the correct number of clusters, as shown in Figure 4b. When using RES distributions, such a high empirical probability of detection cannot be reached, see Figure 4a. As the amount of outliers increases and the performance of the Schwarz BIC starts to deteriorate, the KL divergence remains lower than for the case with the true number of clusters. This can be explained by the degree of freedom of opening up an additional cluster for the outliers, which leads to a deceptive better data fit, that is actually overfitting.
In Figure 5, six exemplary sensitivity curves are shown. The number of clusters is estimated with the Schwarz/skew-Schwarz BIC. As desired for a robust estimator, for this example, the skew-Huber estimator has the lowest KL divergence. It is almost not influenced by the position of the outlier, as indicated by the almost constant KL divergence over the grid. For the Gauß and skew-Gauß estimators, the performance quickly deteriorates as soon as the outlier is positioned outside the clusters. The other robust estimators, i.e. the t, skew-t and Huber show a stable but inferior performance compared to the skew-Huber estimator.
Tables III and IV show the confusion matrices with 0% and 2% of outliers. In the 0% case, all methods have a very high probability of assigning the data points to the correct cluster, with the best result for both Huber (RES and RESK) estimators. For 2% of outlier contamination, both Gaussian estimators break down, whereas the t and Huber based estimators are able to maintain the high performance from the 0% case, again, with the best result for the skew-Huber distribution based estimator.

G. Results Based on Dataset 2: Assessment of Bias and Convergence Speed
In this section, we examine the bias behavior and convergence speed based on Dataset 2. Because this data follows a skew-t distribution with ν = 3, the skew-t, ν = 3 is the MLE for = 0% outliers. In Figure 6, we show the KL divergence over different amounts of outliers. As expected, the best performance can be observed for the skew-t estimator, closely followed by the skew-Huber estimator. As it can be seen from Figure 6b, the KL divergence nearly approaches zero for = 0%, which indicates, in combination with the low variance, that the bias of the algorithms is very small.
To asses the convergence speed of the EM algorithm, the average number of EM iterations is plotted in Figure 7. Two main observations can be made. Firstly, the RESKbased EM algorithm takes much longer to converge than the RES algorithm and, additionally, it has a much more varying stopping time. The anomaly of the skew-Gaussian based EM stopping faster for lager amount of outliers, can be probably  explained by the fact, that the algorithm completely breaks down and stops prematurely at a local optimum.
Lastly, to further investigate the bias behavior and convergence speed, the results of the EM algorithms are shown for different amounts of samples per cluster. No outliers are included in this case. As shown in Figure 8, the larger the amount of samples per cluster, the better the estimation of the true pdf and the lower the variance in the estimation. For a very large number of samples per cluster, the estimated and true pdf are almost identical. This confirms, for the given example, that the approximate ML estimate is near-optimal if the nominal model that is used in the EM algorithm and the data generating model match exactly. For the the RESK algorithms, the convergence speed of the EM algorithms increases for larger sample sizes, as shown in Figure 9.
H. Real Data Simulations 1) Wine Data: The wine quality data set from the UCI Machine Learning Repository (https://archive.ics.uci.edu/ml/ datasets/wine+quality), is composed of 12 different attributes. To distinguish between red and white wine samples, we chose four attributes: volatile acidity, residual sugar, chlorides and total sulfur dioxide, from the data set. The first three attributes are visualized in Figure 2b. The number of clusters, i.e. K = 2 is assumed to be known and, to increase the difficulty, = 3% of uniformly distributed replacement outliers are added in every dimension. From Table V, it is observed that the  modeling of skewness increases the probability of assigning a data point to the correct cluster. Again, the skew-Huber shows the best results.
2) Crabs Data: The crab data set was first used by [62] and can be downloaded at https://doi.org/10.24097/wolfram.70344. data. It is composed of five different length measurements on Leptograpsus crabs, of which all five were used to distinguish between blue/orange and male/female crabs. The number of clusters K = 4 is assumed to be known and no additional outliers are added. Table VI and VII show, that the modeling of skewness, again, increases the overall performance. For this data set, the skew-t distribution based estimator showed best overall performance.

VII. CONCLUSION
We have presented a finite mixture model based on the proposed family of RESK distributions which is capable of integrating robustness and skewness into a single framework. Special attention was given to a newly proposed skew-Huber distribution, which has the attractive property that resulting estimators give maximal weight to data points with a distance smaller than a chosen threshold value while heavily down weighting outliers. The down weighting is asymmetric by taking skewness of the distribution into account. For the proposed family of RESK distributions, we derived an EM algorithm to estimate the cluster parameters and memberships. The performance was evaluated numerically in terms of parameter estimation and clustering performance for skewed and heavy-tailed mixture models, robustness against outliers, and convergence speed in terms of the number of required EM iterations. The proposed methods demonstrated promising performance on all simulated and real-world examples. This showcases their usefulness in simultaneously addressing skewness, outliers and heavy tailed data. The provided framework, therefore, could be useful in a variety of cluster analysis applications where the data may be skewed and heavy-tailed. Future work may consider deriving new robust cluster enumeration criteria and further robust M-estimators to cluster skewed data with outliers. New asymmetric M-estimators can be designed by varying the density generator within the class of RESK distributions according to the user's preference in terms of weighting the data points. Future work may also include deriving a stochastic representation for the skew-Huber distribution or other M-estimators. Based on this stochastic representation, the corresponding approximate ML estimator could be replaced by the exact ML estimator.

APPENDIX A DERIVATION OF M-STEP ESTIMATES
Using the matrix calculus rules from [56], [58], F is a differentiable m × p matrix function of a n × q matrix X. Then, the Jacobian matrix of F at X is a mp × nq matrix DF (X) = ∂vec(F (X)) ∂(vec(X)) .
Starting with the differential of the complete data loglikelihood regarding ξ m , we define F as a 1×1 scalar function of the r × 1 vector ξ m . Hence, the resulting Jacobian matrix is just the gradient and is of size 1 × r. Setting F equal to (34), i.e., 2 Ω −1 Both differentials are calculated separately as with the Sherman-Morrison formula [63] A + bc and dF (κ nm ) with Using the scalar values the differentials can be combined to the gradient as follows Setting Equation (55) to zero and solving for ξ m leads to the ML estimate. As discussed in Section IV, finding an exact solution, is not feasible for general RESK distributions. Therefore, the values e 0,nm and e 1,nm , are replaced by their approximations, i.e.ẽ 0,nm ,ẽ 1,nm . These approximations are constants with respect to the unknown parameters θ m . They are introduced by approximating t nm , τ m , η nm , and κ nm with their corresponding estimates from the previous iteration of Algorithm 1. The resulting approximate ML estimate is the obtained by setting the approximate gradient to zero: Now, F is defined as a 1×1 scalar function of the r ×1 vector λ m . Hence, the resulting Jacobian matrix is just the gradient and is of size 1 × r. Setting F equal to (29) and applying the differential gives Again, the occurring differentials are calculated separately. Firstly, with dvec λ m λ m = ((λ m ⊗ I r ) + (I r ⊗ λ m )) dλ m . (59) with vectorization dvec (g (t nm )) and thirdly with vectorization Combining the differentials with e 2,nm =τ 2 m + ψ (t nm ) 2η 2 nm − Ψ (κ nm ) η nm τ m 2ψ(t nm ) Analogously to Eq. (55), the gradient in Eq. (65) is approximated by replacing e 1,nm and e 2,nm by their approximations, i.e.ẽ 1,nm ,ẽ 2,nm . Then, the approximate ML estimate follows by setting the approximate gradient to zero and solving for λ m :λ Lastly, F is defined as a 1 × 1 scalar function of the r × r matrix S m . Hence, the resulting Jacobian matrix is just the gradient and is of size 1 × r 2 . Setting F equal to (29) and applying the differential dF (S m ) =  dg (t nm ) F (κ nm ) again, leads to three differentials that are calculated separately. Starting with Here, we used the following properties of the commutation and duplication matrix K n,n = K n , K n D n = D n (73) where A ∈ R m×n , B ∈ R p×q , K m,n ∈ R mn×mn and b ∈ R p×1 . Combining the results leads to the gradient DF (S m ) = Analogously to Eq. (55) and Eq. (65), the gradient in Eq. (75) is approximated by replacing e 0,nm , e 1,nm and e 2,nm by their approximations, i.e.ẽ 0,nm ,ẽ 1,nm ,ẽ 2,nm . Then, the approximate ML estimate follows by setting the approximate gradient to zero and solving for S m : and, after the elimination of λ, we find ACKNOWLEDGMENT Christian A. Schroth is supported by the DFG Project Number 431431951. The work of Michael Muma has been funded by the LOEWE initiative (Hesse, Germany) within the emergenCITY centre and is supported by the 'Athene Young Investigator Programme' of Technische Universität Darmstadt, Hesse, Germany.