## 3 Material and Methods
### MRCD Algorithm  
The Minimum Covariance Determinant (MCD) based method is one of the most extensively employed robust multivariate location and scatter estimators. This method is designed to determine \(\mathbf{H}_{MCD}\), defined as a subset with the smallest sample covariance determinant. The MCD estimates for the mean vector and covariance matrix correspond to the mean vector and covariance matrix of \(\mathbf{H}_{MCD}\). Define \(h\) as the subset size, where \(\frac{n}{2}\leq h<n\) and \(h\geq p\) must be fulfilled; otherwise, the covariance matrix derived from the MCD method would become singular. The MCD algorithm calculates every subset  
Figure 1: Bootstrap Control Limit Algorithm. possible, as many as \( \begin{pmatrix}n\\ h\end{pmatrix} \) possible combinations, in order to obtain \(\mathbf{H}_{MCD}\). Therefore, this method is time-consuming and not suitable for estimating large datasets. Minimum Regularized Covariance Determinant (MRCD) estimators are proposed [23] as the extension of MCD. The MRCD is a robust estimator that uses various combinations of target matrix and regularization weight determined through data-driven procedures. The application of the MRCD estimator is robust and well-suited for handling outliers in datasets with high dimensions. Let \(\mathbf{X}=\left\{\mathbf{x}_{1},\,\mathbf{x}_{2},\,\ldots,\,\mathbf{x}_{i},\ldots,\,\mathbf{x}_{n }\right\}^{\prime}\), where \(\mathbf{x}_{i}=\left\{x_{i1},\,x_{i2},\,\ldots,\,x_{ip}\right\}^{\prime}\) from a \(p\)-variate observations. First, the data must be standardized using the Qn estimator [41] and these values must be put in a diagonal matrix \(\mathbf{D}_{X}\). The median of each variable also needs to be computed and put in a location vector \(\mathbf{v}_{X}\). The standardized observations are then stated on \(\mathbf{U}\) that constructed under a set of \(\mathbf{u}_{i}\) as follows:  
\[\mathbf{u}_{i}=\mathbf{D}_{X}^{-1}(\mathbf{x}_{i}-\mathbf{v}_{X}) \tag{3}\]  
The next step is defining the target matrix (\(\mathbf{\Gamma}\)) and scalar regularization parameter \((\rho)\). \(\mathbf{\Gamma}\) is a \(p\times p\) diagonal matrix that consists of estimated univariate scales, while \(\rho\) is a weighted parameter that can be obtained by a data-driven approach that satisfies \(0\leq\rho\leq 1\). Then, define \(h\times p\) covariance matrix \(\mathbf{K}(H)\) of \(h\)-subset of \(H\) on the standardized data \(\mathbf{U}\) as follows:  
\[\mathbf{K}(H)=\rho\;\mathbf{\Gamma}+(1-\rho)c_{\gamma}\mathbf{S}_{\mathbf{U}}(H) \tag{4}\]  
where \(\mathbf{S}_{\mathbf{U}}(H)\) is the covariance matrix of the \(h\)-subset for \(\mathbf{U}\) and \(c_{\gamma}\) is the consistency factor [42]. Mathematical operation (7) can be done by a spectral decomposition \(\mathbf{\Gamma}=\mathbf{Q}\mathbf{\lambda}\mathbf{Q}^{\prime}\) where \(\mathbf{\lambda}\) represents the diagonal matrix encompassing the eigenvalues of \(\mathbf{\Gamma}\), while \(\mathbf{Q}\) denotes the orthogonal matrix comprising the corresponding eigenvectors. The previous equation can be rewritten as follows:  
\[\mathbf{K}(H)=\mathbf{Q}\mathbf{\lambda}^{1/2}[\rho\;\mathbf{I}+(1-\rho)c_{\alpha}\mathbf{S}_\mathbf{{W}}(H)]\mathbf{\lambda}^{1/2}\mathbf{Q}^{\prime} \tag{5}\]  
where \(\mathbf{W}\) contains the transformation of standardized observations \(\mathbf{w}_{i}=\mathbf{\lambda}^{1/2}\mathbf{Q}^{\prime}\mathbf{u}_{i}\). Consequently, it follows \(\mathbf{S}_{\mathbf{W}}(H)=\mathbf{\lambda}^{-1/2}\mathbf{Q}^{\prime}\mathbf{S}_{\mathbf{U}}(H)\mathbf{Q }\mathbf{\lambda}^{-1/2}\).