New Highly Efficient High-Breakdown Estimator of Multivariate Scatter and Location for Elliptical Distributions

The elliptical distribution is a general class of multivariate probability distributions encompassing many familiar subclasses such as the symmetric Gaussian, t-, Cauchy, Laplace, hyperbolic, variance gamma, and normal inverse Gaussian distributions. Table 1 summarizes the most common elliptical distributions (Fang et al., 1990, p. 69; Deng and Yao, 2018).

Symmetric elliptical distributions are defined as being a function of the squared Mahalanobis distance,¹¹1Some texts define the Mahalanobis distance with the mean and covariance, but this more restrictive definition excludes thick-tailed distributions where these do not exist, such as Cauchy distributions. $d\left(\boldsymbol{x},\boldsymbol{\mu},\boldsymbol{\Sigma}\right)=\left(\boldsymbol{x}-\boldsymbol{\mu}\right)^{\operatorname{T}}\boldsymbol{\Sigma}^{-1}\left(\boldsymbol{x}-\boldsymbol{\mu}\right)$, where $\boldsymbol{x}\in\mathbb{R}^{p}$, the location $\boldsymbol{\mu}\in\mathbb{R}^{p}$, and the $p\times p$ positive definite symmetric ($\operatorname{PDS}(p)$) scatter $\boldsymbol{\Sigma}\in\operatorname{PDS}(p)$. When the PDF is defined, it has the form $f_{X}\left(\boldsymbol{x}\right)=\alpha_{p}|\boldsymbol{\Sigma}|^{-1/2}\phi\left(d\left(\boldsymbol{x},\boldsymbol{\mu},\boldsymbol{\Sigma}\right)\right),$ for some generating function $\phi(d)$, and where $\alpha_{p}$ is a constant that ensures $f_{X}\left(\boldsymbol{x}\right)$ integrates to one. Table 1 lists common generating functions. When the covariance exists, it is proportional to the scatter matrix, $\boldsymbol{\Sigma}$. The corresponding shape matrix is commonly defined as

The PDF of $d\left(\boldsymbol{x},\boldsymbol{\mu},\boldsymbol{\Sigma}\right)$ is given by (Kelker, 1970)

where $\beta_{p}=\alpha_{p}\pi^{p/2}/\Gamma(p/2)$. Hereafter, all densities, $f(d),$ refer to the density of $d\left(\boldsymbol{x},\boldsymbol{\mu},\boldsymbol{\Sigma}\right)$ in (2), so the subscript $D$ will be omitted. It is also generally assumed that $p>2.$

Given a set of $n$ $p$-dimensional samples, $\{\boldsymbol{x}_{1},...,\boldsymbol{x}_{n}\}$, S-estimators of location and shape are defined as (Maronna et al., 2006, Sec. 6.4.2)

for some scalar rho function, $\rho(t)$. A proper S-estimator rho function should be a continuously differentiable, nondecreasing function in $t\geq 0$ with $\rho(0)=0$, and where there is a point $c$ such that $\rho(t)=\rho(\infty)$ for $t\geq c$. For simplicity, and without loss of generality, the rho functions will be normalized so $\rho(\infty)=1$. The parameter $b$ is a scalar that affects the efficiency (see Section 4) and robustness of the estimator. The purpose of S-estimators is to achieve high robustness, so they are usually configured with $b=1/2-(p+1)/(2n),$ which achieves the maximum theoretical breakdown point that any affine equivariant estimator may have (see Section 5.1). To understand the derivation of the proposed estimator in the next section, note that $\widehat{\sigma}$ in the constraint is an M-estimator of the scale of $d\left(\boldsymbol{\mu},\boldsymbol{\Omega}\right).$ Local solutions of (3) can be found iteratively using the weighted sums $\sum_{i=1}^{n}w\left(d_{i}/\widehat{\sigma}\right)\left(\boldsymbol{x}_{i}-\widehat{\boldsymbol{\mu}}\right)=\boldsymbol{0}$ and $\sum_{i=1}^{n}w\left(d_{i}/\widehat{\sigma}\right)\left(\boldsymbol{x}_{i}-\widehat{\boldsymbol{\mu}}\right)\left(\boldsymbol{x}_{i}-\widehat{\boldsymbol{\mu}}\right)^{\operatorname{T}}\propto\widehat{\boldsymbol{\Omega}},$ where the weight function $w(t)=\rho^{\prime}(t),$ and where $\widehat{\boldsymbol{\Omega}}$ is re-normalized with each iteration. For the empirical results in this paper, the estimators will all be solved using this weighted-sum algorithm.

To estimate the scatter metrix, a separate estimator of $|\boldsymbol{\Sigma}|^{1/p}$ can then be used to scale $\widehat{\boldsymbol{\Omega}}$ using (1). Maronna et al. (2006, p. 186) discussed a simple estimator to scale $\widehat{\boldsymbol{\Omega}}$ to $\widehat{\boldsymbol{\Sigma}}$. When $\boldsymbol{x}$ is normally distributed, $d$ has a chi-squared distribution with $p$ degrees of freedom. Therefore, they suggested using $\widehat{\boldsymbol{\Sigma}}=\operatorname{Median}\left\{d\left(x_{1},\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Omega}}\right),\dots,d\left(x_{n},\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Omega}}\right)\right\}\left(\chi^{2}_{p}(0.5)\right)^{-1}\widehat{\boldsymbol{\Omega}}$, where $\chi^{2}_{p}(0.5)$ is the 50^th percentile of the chi-squared distribution. For the general case of elliptical distributions, we propose extending this to

where $F(d)$ is the distribution function corresponding to (2), and therefore $F^{-1}(0.5)$ is the 50^th percentile of the distribution.

For the location and shape matrices, the S-estimator formulation in (3) is equivalent to the alternative one given by

which requires that $\sigma$ be defined such that $b=\operatorname{E}\left[\rho\left(d\left(\boldsymbol{x};\boldsymbol{\mu},\boldsymbol{\Sigma}\right)/\sigma\right)\right]$ for a consistent estimator of $\boldsymbol{\Sigma}$ at an assumed elliptical distribution (Rocke, 1996). While the first formulation is better for understanding the derivation of the proposed S-q estimator, this second formulation is better for defining and understanding its properties (for example, see Lopuhaä, 1989). The scale parameters in the two formulations are related asymptotically by $\sigma=|\boldsymbol{\Sigma}|^{-1/p}\operatorname{E}\left[\widehat{\sigma}\right],$ at the assumed distribution.

The two most common multivariate S-estimators are the bisquare and Rocke (Maronna et al., 2019, sec. 6.4.2, 6.4.4). The S-bisquare is given by $\rho_{\text{bisq}}(t)=\min\left\{1,1-(1-t)^{3}\right\}$ and $w_{\text{bisq}}(t)=3(1-t)^{2}\operatorname{I}(t\leq 1),$ which does not have a tuning parameter to control efficiency and robustness. The S-Rocke is given by

where the parameter $\gamma\in(0,1]$ tunes the estimator’s efficiency and robustness. The Rocke’s maximum efficiency is generally limited at $\gamma=1$, which is extremely restricting for small $p$. Both $\rho_{\text{bisq}}(t)$ and $\rho_{\gamma}(t)$ are generic functions that do not depend on the underlying distribution. In the following section, an alternative S-estimator is defined that accounts for the underlying distribution and that generally has better performance across the most common elliptical distributions. It also does not have the same inherent restrictions for small $p$ as $\rho_{\gamma}(t).$

The rho function corresponding to the maximum likelihood estimator, $\widehat{\sigma},$ of the scale of $d\left(\boldsymbol{\mu},\boldsymbol{\Omega}\right),$ or equivalently $d\left(\boldsymbol{\mu},\boldsymbol{\Sigma}\right),$ is $\rho_{\text{mle}}(t)=-t\,f^{\prime}\left(t\right)/f\left(t\right).$ We propose weighting this by the power transform of the density, $\tilde{\rho}_{q}\left(t\right)=f(t)^{1-q}\rho_{\text{mle}}(t),$ where the scalar $q\leq 1$ controls the estimator robustness, with $q=1$ corresponding to the maximum likelihood function, and with decreasing $q$ increasing the estimator robustness. In most cases, this rho function is not monotone, as required by S-estimators, so it is denoted with a tilde. This rho function is equivalent to the M-Lq and other M-estimators proposed, for example, by Windham (1995); Basu et al. (1998); Choi and Hall (2000); and Ferrari and Yang (2010). However, in this particular case of estimating the scale of the squared Mahalanobis distance of an elliptically distributed random vector, the density and rho function do not need to be regenerated with each numerical iteration, $i,$ based on the estimates $\widehat{\boldsymbol{\mu}}^{(i)}$ and $\widehat{\boldsymbol{\Omega}}^{(i)}.$ Substituting the PDF from (2),

where $s_{p}=p/2-1$ and $s_{q}=1-q$. Taking the derivative of $\tilde{\rho}_{q}(t)$, the corresponding weight function is given by

For simplicity, when $q<1$, the scalar $\beta_{p}$ can be dropped from the calculation of $\tilde{\rho}_{q}(t)$ and $\tilde{w}_{q}(t)$ in (5) and (6). When $\phi(t)$ is only positive over a finite domain (e.g. Pearson Type II distribution), then we define $\tilde{\rho}_{q}(t)$ and $\tilde{w}_{q}(t)$ to be zero outside this domain.

For the common elliptical distributions listed in Table 1, $\tilde{\rho_{q}}(t)$ is monotone in its central region between its global extrema when using appropriate values for $q$ (defined below). The first extremum is the minimum, which we label point $a$, and the second is the maximum, labeled $c$. The distance between $a$ and $c$ varies monotonically with respect to $q$. We use this to define a tunable, double-hard-rejection S-estimator rho function. The value of $\tilde{\rho}_{q}(t)$ is held constant between zero and $a$ at value $\tilde{\rho}_{q}(a),$ which hard rejects inliers, and the value of $\tilde{\rho}_{q}(t)$ is held constant above $c$ at value $\tilde{\rho}_{q}(c),$ which hard rejects outliers. The resulting monotonic function is then scaled and shifted so it ranges from zero to one. This defines the S-q estimator.

For most distributions, $\lim_{q\rightarrow 1}c=\infty,$ or at $q=1$, $\tilde{\rho}_{q}(t)$ is not bounded. Therefore, we do not scale or shift $\tilde{\rho}_{q}(t)$ in this case, and $\rho_{q}(t)$ is not a proper S-estimator rho function. However, when $q=1$ and $b=1$, the MLE of the scale of $d$ is obtained. The S-q weight function is the derivative of $\rho_{q}(t)$ and is given by

Table 2 lists expressions for the inlier rejection point, $a,$ and the outlier rejection point, $c$, for the common elliptical distributions in Table 1. For most of these distributions, the equation $\tilde{w}_{q}(t)=0$ is quadratic, which provides a closed-form solution for the values of $a$ and $c$.

The asymptotic rejection probability (ARP) is defined as $Pr(d/\widehat{\sigma}\geq c)$ (Rocke, 1996). Table 2 can be used to determine $q$ from a desired ARP using $F^{-1}(ARP)$. However, since $w_{q}(t)$ is very tapered (i.e. applying little weight to values just below $c$), practitioners may choose alternative approaches to tuning that allow for higher estimator efficiencies. For example, the approach used in this paper as well as in Maronna and Yohai (2017) is to tune the estimators to a desired expected efficiency, which is defined in the next section.

The general definition in (7) specifies that $q\leq 1$. In a few particular cases, however, there are some minor restrictions on $q$ (when $q<1$) in order to ensure that $a$ and $c$ are in the support of $f(d)$. Table 3 lists these restrictions.

Figure 1 illustrates examples of the S-q functions $\tilde{\rho}_{q}(t)$, $\rho_{q}(t)$, and $w_{q}(t)$ for the five-dimensional Gaussian (S-q Type II) and Laplace (S-q Type I) distributions and for various values of $q$. As $q$ is decreased, the region of positive weights (area between points $a$ and $c$) narrows, corresponding to increased robustness. The PDF is also plotted, illustrating how $w_{q}(t)$ roughly follows $f(t)$ in the central region.

Figure 2 compares the S-q asymptotic weights with those of the MM-SHR, S-Rocke, S-bisquare, and maximum likelihood estimators, and with the corresponding PDF. The underlying model is a 10-dimensional standard Gaussian distribution. The MM-SHR and S-q estimators have been tuned to 80% asymptotic efficiency relative to the MLE. The S-Rocke estimator is tuned to its maximum efficiency, which is 77% in this instance. The estimators have been set to the maximum breakdown point, with $b=1/2,$ which results in the shifts of the peaks of the weight curves relative to the PDF.

From the figure, it is clear that the Gaussian MLE (i.e. sample estimator) gives uniform weight to all samples, no matter how improbable. The S-Rocke has a quadratic weight function, which is a hard cutoff that cannot capture the tails of $f(d)$. The SHR weight function is cubic, and its shape better captures the shape of the right-half of the PDF. However, the SHR function is designed to approximate a step function, which is poorly suited for many distributions (c.f. $w_{shr}(t)$ in Figure 2 with the Laplace $f(d)$ in Figure 1). Only the S-q weight function follows the general shape of the PDF—giving less weight to less probable observations.

For the asymptotic distribution of the S-q estimate, we continue to use the alternative S-estimator formulation given by (4). Lopuhaä (1997) derived the distribution of S-estimators with assumptions appropriate for the S-q estimator, that is

Here, we use the following notation. The matrix $\boldsymbol{I}_{p^{2}}$ is the $p^{2}\times p^{2}$ identity matrix, $\boldsymbol{K}_{p^{2}}$ is the $p^{2}\times p^{2}$ commutation matrix, $\otimes$ is the Kronecker product operator, and the operator $\operatorname{vec}\left(\boldsymbol{\Sigma}\right)$ stacks the columns if $\boldsymbol{\Sigma}$ into a column vector.

Frahm (2009) derived the asymptotic distribution of shape matrix estimates for affine equivariant estimators. This enables us to state the asymptotic distribution of the S-q shape estimate, which is applicable using either S-estimator formulation, (3) or (4).

The asymptotic efficiency of an estimator, at an assumed distribution, is defined as the ratio of the asymptotic variance of the maximum likelihood estimate to the variance of the estimator under consideration. For multivariate estimation, this definition of efficiency is of large dimension—${p\times p}$ for location and ${p^{2}\times p^{2}}$ for shape and scatter. However, for affine equivariant estimation of location and scatter of elliptical distributions, the covariance of the estimate depends only on a scalar. Specifically, (9), (10), and (11) are general, with only the scalars $\omega_{1}/\omega_{2}^{2}$ (Bilodeau and Brenner, 1999), and $\zeta_{1}$ and $\zeta_{2}$ (Tyler, 1982) depending on the estimator and the generating function $\phi(d)$. Therefore, the asymptotic efficiency of the estimate $\widehat{\boldsymbol{\mu}}$ can be alternatively defined as

and the asymptotic efficiency of the estimate $\widehat{\boldsymbol{\Omega}}$ can alternatively be defined as

It is common to define asymptotic efficiency this way (for example, see Tyler, 1983; Frahm, 2009).

Comparing the S-q estimator’s efficiency to another estimator can likewise be achieved analytically using, for example, $\zeta_{1,\gamma}/\zeta_{1,q}$ for the S-Rocke estimator, which when the quotient is greater than one, indicates that the S-q has higher asymptotic efficiency than the S-Rocke estimator. For other S-estimators, the asymptotic distribution parameters $\omega_{1},$ $\omega_{2},$ and $\zeta_{1}$ are calculated the same as in Theorems 4 and 5 but using their respective weight functions. MM-estimators have the same asymptotic variance and influence function as S-estimators (Rousseeuw and Hubert, 2013). For MM-estimators, however, $\sigma$ is effectively the tuning parameter, and it can be set accordingly.

In general, finite-sample performance measures are difficult to derive analytically. Instead, it is common to characterize finite-sample performance by empirically characterizing the behavior of metrics derived from the Kullback-Leibler divergence between the estimated and true distribution (for example, see Huang et al., 2006; Ferrari and Yang, 2010). For t-distributions, which includes the Gaussian distribution, the Kullback-Leibler divergence between $t_{\nu}\left(\boldsymbol{\mu},\boldsymbol{\Sigma}\right)$ and $t_{\nu}\left(\widehat{\boldsymbol{\mu}},\widehat{\boldsymbol{\Sigma}}\right)$ is given by (Abusev, 2015)

Following Maronna and Yohai (2017), we then define the joint location and scatter finite-sample relative efficiency as

where $\widehat{\boldsymbol{\mu}}_{\text{mle}}$ and $\widehat{\boldsymbol{\Sigma}}_{\text{mle}}$ are the location and scatter matrices corresponding to the maximum likelihood estimate, and where the expectation is calculated empirically using the sample mean over $m$ Monte Carlo trials. The location and the scatter finite-sample relative efficiencies are then respectively defined as $\operatorname{eff}_{n}\left(\widehat{\boldsymbol{\mu}},\boldsymbol{\Sigma};\widehat{\boldsymbol{\mu}}_{\text{mle}},\boldsymbol{\Sigma}\right)$ and $\operatorname{eff}_{n}\left(\boldsymbol{\mu},\widehat{\boldsymbol{\Sigma}};\boldsymbol{\mu},\widehat{\boldsymbol{\Sigma}}_{\text{mle}}\right).$ Likewise, we define the shape matrix finite-sample relative efficiency as

Any estimator must provide a good estimate in the absence of contamination and when tuned to its maximum efficiency. This section compares the maximum achievable efficiencies of the S-q, S-Rocke, and MM-SHR estimators when set to their maximum breakdown point. The results below cover large swaths of the most common elliptical families in Table 1 for a moderate dimension of $p=20.$ These swaths were specifically chosen to cover everyday distributions: Gaussian, Cauchy, Laplace, hyperbolic, and normal inverse Gaussian distributions.

Robust scatter matrix estimation is generally “more difficult” than the estimation of location (Maronna et al., 2019), and as Maronna and Yohai (2017) demonstrated, divergence and efficiency metrics for scatter matrix estimators are generally much worse than for the corresponding estimators of location. Likewise, due to the high dimensionality of the estimate, the underlying shape matrix is the most difficult part of estimating the scatter matrix. Additionally, many practical applications such as multivariate regression, principal components analysis, linear discriminant analysis, and canonical correlation analysis only require the shape matrix, and not the full scatter or covariance matrices (Frahm, 2009). Therefore, unless otherwise noted, the performance results in this paper are for the shape matrix, with metrics given by (12), (13), and $D\left(\boldsymbol{\mu},\boldsymbol{\Omega};\boldsymbol{\mu},\widehat{\boldsymbol{\Omega}}\right).$

The maximum efficiencies of the S-q and S-Rocke generally occur when their parameters $q$ and $\gamma$ are set to one—although the maximum breakdown point of the S-q is only achieved when $q<1$. However, the maximum efficiency of the MM-SHR must be determined by a search as depicted in Figure 3, which plots, as an example, asymptotic efficiency versus tuning parameter for the estimators for the 20-dimensional Cauchy distribution. At the limit, as the MM-SHR parameter is increased toward infinity, all samples receive equal weight, which is the MLE for the Gaussian distribution, but not for distributions such as the Cauchy. In general, for each tunable estimator, its efficiency decreases while its robustness increases as its parameter is decreased. At the lower limit of its parameter, its weight function is a delta function that may reject all the samples and may result in zero efficiency. At this point, the robustness is high, but the weighted-sum solution depends entirely on the initial estimates $\widehat{\boldsymbol{\mu}}^{(0)}$ and $\widehat{\boldsymbol{\Omega}}^{(0)}.$

It should be noted that although generally of high efficiency, the S-q estimate at its limit with $q=1$ is not necessarily the maximum likelihood estimate for location and scatter. The MLE weight function for location and scatter is given by (Tyler, 1982)

whereas at $q=1,$ (8) gives

An example family that satisfies this theorem is the Kotz type with parameter $N=0.$ Note, however, that $\lim_{q\rightarrow 1}\tilde{\rho}_{q}(c)=\infty,$ so high breakdown cannot be achieved simultaneously. This is illustrated in Figure 4, which provides the estimators’ maximum achievable asymptotic shape efficiencies for the Kotz type distribution with parameters $N=0$ and $r=1/2$ as a function of parameter $s$. In this example, the S-q efficiency is plotted for its maximum absolute efficiency with $q=1$ and for it approximate maximum high-breakdown efficiency with $q=0.99.$ As seen in the figure, the high cost of high-breakdown is particularly acute for large $s.$

The remainder of this paper will focus on maximum efficiency at the maximum breakdown point. Figure 4 also provides the S-Rocke and MM-SHR estimator’s maximum efficiencies at their maximum breakdown points. The MM-SHR efficiency peaks at $s=1,$ which is expected since this is the Gaussian distribution, and the S-Rocke efficiency peaks just above this point. Their efficiencies fall off precipitously for larger and smaller values of $s.$ The efficiency of the S-q, conversely, increases toward unity for smaller $s.$

The estimators’ maximum achievable asymptotic efficiencies for the t-distribution as a function of the distribution parameter, $\nu$, are plotted on the left of Figure 5. When $\nu=1$, the t-distribution corresponds to a Cauchy distribution, and when $\nu\rightarrow\infty,$ it corresponds to the Gaussian distribution. The S-q estimator offers the highest efficiency of the three estimators for thicker tails.

The maximum achievable small-sample relative efficiencies using $n=3p$ are plotted on the right of Figure 5. The initial estimates were made using the Peña and Prieto (2007) kurtosis plus specific directions (KSD) estimator as recommended and provided by Maronna and Yohai (2017). Comparing these finite-sample results with the asymptotic ones on the left, it is seen that the relative results are similar. This general similarity implies that the relative performance of the asymptotic efficiencies can often be a good surrogate for the relative performance of the finite-sample efficiencies when there is no closed-form expression for the divergence in (13).

The estimators’ maximum achievable asymptotic efficiencies for the variance gamma distribution with $\psi=2$ are plotted as a function of parameter $\lambda$ on the left of Figure 6. The plots highlight the Laplace $(\lambda=1)$ and multivariate hyperbolic $(\lambda=(p+1)/2)$ distributions. The S-q exhibits good performance for the hyperbolic and remarkably good performance for the Laplace.

The estimators’ maximum achievable asymptotic efficiencies for the generalized hyperbolic distribution with $\psi=2$ and $\chi=1$ are plotted as a function of parameter $\lambda$ on the right of Figure 6. The plots highlight the normal inverse Gaussian $(\lambda=-1/2)$ and hyperbolic $(\lambda=(p+1)/2)$ distributions. The S-q again exhibits good performance for the hyperbolic, and it exhibits remarkably good performance for the normal inverse Gaussian.

The finite-sample breakdown point of a multivariate estimator of location or scatter is defined as the fraction of the samples, $\epsilon n,$ that can be set to either drive $\|\widehat{\boldsymbol{\mu}}\|=\infty$ or drive an eigenvalue of $\widehat{\boldsymbol{\Sigma}}$ to either zero or infinity. Unlike multivariate M-estimators, which only achieve an asymptotic breakdown point of $(p+1)^{-1}$ (Maronna, 1976), S-estimators are able to achieve the maximum possible finite-sample breakdown point that any affine equivariant estimator may have (Davies, 1987, Th. 6). For the following theorem, the term samples in general position means that no more than $p$ samples are contained in any hyperplane of dimension less than $p$.

The influence function (IF) of an estimator characterizes its sensitivity to an infinitesimal point contamination at $\boldsymbol{z}\in\mathbb{R}^{p},$ standardized by the mass of the contamination, $\epsilon.$ The influence function for estimator $\boldsymbol{T},$ at the nominal distribution $F$, is defined as

where $\epsilon$ is the proportion of samples that are a point-mass, $\Delta_{\boldsymbol{z}},$ located at $\boldsymbol{z}.$

By definition of S-estimators with normalized rho function, the magnitude of first term of (17) is clearly bounded to no more than $\lambda_{2}^{-1}\boldsymbol{\Sigma}.$ Therefore, to compare the influence functions of the S-q, S-Rocke, and MM-SHR estimators, we focus on the second term. From this term, define $\alpha_{\boldsymbol{\Sigma}}\left(d_{z}\right)=\lambda_{1}^{-1}p(p+2)(d_{z}/\sigma)w(d_{z}/\sigma)$ for each estimator. Figure 7 plots $\alpha_{\boldsymbol{\Sigma}}(d_{z})$ at the 10-dimensional Gaussian distribution for the estimators as depicted in Figure 2.

By definition, all highly-robust estimators have bounded influence functions, and for the three estimators considered here, their influence functions are continuous. This means that small amounts of contamination have limited effects on their estimates. The gross-error sensitivity of an estimator is the maximum of $\boldsymbol{\operatorname{IF}}\left(\boldsymbol{z}\right),$ and in this example, the S-q demonstrates a lower gross-error sensitivity than the S-Rocke and MM-SHR estimators. By its definition, the MM-SHR has a inlier rejection point of zero, meaning inliers can negatively influence its estimates. However, proper Type II S-q functions have positive inlier rejection points, which provide robustness against inliers.

Relative to the S-Rocke and MM-SHR estimators, the S-q often has larger outlier rejection points. This is the cost of its generally higher efficiency and ability to reject inliers. However, due to its continuity, the influence near this point is still greatly attenuated.

To empirically compare the finite-sample robustness of the estimators, we employed the simulation method used by Maronna and Yohai (2017) and plot the shape matrix divergence, $\operatorname{D}\left(\boldsymbol{\mu},\boldsymbol{\Omega};\boldsymbol{\mu},\widehat{\boldsymbol{\Omega}}\right)$, versus shift contamination value $k$. For a contamination proportion $\epsilon$, the first element of each of the $\lfloor\epsilon n\rfloor$ contaminated samples was replaced with the value $k$, that is $x_{1}=k$. The initial estimates of the weighted algorithm were determined with the KSD estimator. Figure 8 provides divergence plots for normally distributed data with $\epsilon=10\%$ contamination, for dimensions $p=5$ and $p=20$, and for sample sizes $n=5p$ and $n=100p$. For the cases where $p=20$, the estimators were tuned to 90% uncontaminated relative efficiency. When $p=5$, the S-Rocke has poor maximum efficiency, so the estimators were tuned to match the maximum S-Rocke efficiency.