This vignette illustrates the usage of the package fitHeavyTail to estimate the mean vector and covariance matrix of heavy-tailed multivariate distributions such as the angular Gaussian, Cauchy, or Student’s \(t\) distribution. The results are compared against existing benchmark functions from different packages.


The package can be installed from CRAN or GitHub:

# install stable version from CRAN

# install development version from GitHub

To get help:

help(package = "fitHeavyTail")

To cite fitHeavyTail in publications:


Quick Start

To illustrate the simple usage of the package fitHeavyTail, let’s start by generating some multivariate data under a Student’s \(t\) distribution with significant heavy tails (degrees of freedom \(\nu=4\)):

library(mvtnorm)  # package for multivariate t distribution
N <- 10   # number of variables
T <- 80   # number of observations
nu <- 4   # degrees of freedom for heavy tails

mu <- rep(0, N)
U <- t(rmvnorm(n = round(0.3*N), sigma = 0.1*diag(N)))
Sigma_cov <- U %*% t(U) + diag(N)  # covariance matrix with factor model structure
Sigma_scatter <- (nu-2)/nu * Sigma_cov
X <- rmvt(n = T, delta = mu, sigma = Sigma_scatter, df = nu)  # generate data

We can first estimate the mean vector and covariance matrix via the traditional sample estimates (i.e., sample mean and sample covariance matrix):

mu_sm     <- colMeans(X)
Sigma_scm <- cov(X)

Then we can compute the robust estimates via the package fitHeavyTail:

fitted <- fit_mvt(X)

We can now compute the estimation errors and see the significant improvement:

sum((mu_sm     - mu)^2)
#> [1] 0.2857323
sum((fitted$mu - mu)^2)
#> [1] 0.1487845

sum((Sigma_scm  - Sigma_cov)^2)
#> [1] 5.861138
sum((fitted$cov - Sigma_cov)^2)
#> [1] 3.031499

To get a visual idea of the robustness, we can plot the shapes of the covariance matrices (true and estimated ones) on two dimensions. Observe how the heavy-tailed estimation follows the true one more closely than the sample covariance matrix:

Numerical Comparison with Existing Packages

In the following, we generate multivariate heavy-tailed Student’s \(t\) distributed data and compare the performance of many different existing packages via 100 Monte Carlo simulations in terms of estimation accurary, measured by the mean squared error (MSE) and CPU time.

The following plot gives a nice overall perspective of the MSE vs. CPU time tradeoff of the different methods (note the ellipse at the bottom left that embraces the best four methods: fitHeavyTail::fit_Tyler(), fitHeavyTail::fit_Cauchy(), fitHeavyTail::fit_mvt(), and fitHeavyTail::fit_mvt() with fixed nu = 6):

From the numerical results we can draw several observations:

Concluding, the top choices seem to be (in order):

  1. fitHeavyTail::fit_mvt() (either without fixing nu or with nu = 6),
  2. fitHeavyTail::fit_Cauchy(),
  3. fitHeavyTail::fit_Tyler(), and
  4. MASS::cov.trob() (with the advantage of being preinstalled with base R, but with a worse estimation error).

The overall winner is fitHeavyTail::fit_mvt() by a big margin.

Extension to Skewed Distributions

The empirical distribution of daily returns of some financial variables, such as exchange rates, equity prices, and interest rates, is often skewed. There are several different formulations of multivariate skewed \(t\) distributions appearing in the literature (Lee and McLachlan, 2014) (Aas and Haff, 2006). The package now supports the multivariate generalized hyperbolic (GH) skewed \(t\) distribution and provides a method to estimate the parameters of such distribution. It is implemented in the function fitHeavyTail::fit_mvst(). Below is a simple example to illustrate its usage:

# parameter setting for GH Skewed t distribution
N <- 5
T <- 200
nu <- 6
mu <- rnorm(N)
scatter <- diag(N)
gamma <- rnorm(N)

# generate GH Skew t data via hierarchical structure
taus <- rgamma(n = T, shape = nu/2, rate = nu/2)
X <- matrix(data = mu, nrow = T, ncol = N, byrow = TRUE) +
     matrix(data = gamma, nrow = T, ncol = N, byrow = TRUE) / taus +
     mvtnorm::rmvnorm(n = T, mean = rep(0, N), sigma = scatter) / sqrt(taus)

# fit GH Skew t model
fitted <- fit_mvst(X)


In essence, all the algorithms are based on the maximum likelihood estimation (MLE) of some assumed distribution given the observed data. The difficulty comes from the fact that the optimal solution to such MLE formulations becomes too involved in the form of a fixed-point equation and the framework of Majorization-Minimization (MM) algorithms (Sun et al., 2017) becomes key to derive efficient algorithms.

In some cases, the probability distribution function becomes too complicated to manage directly (like the multivariate Student’s \(t\) distribution) and it is necessary to resort to a hierarchical distribution that involves some latent variables. In order to deal with such hidden variables, one has to resort to the Expectation-Maximization (EM) algorithm, which interestingly is an instance of the MM algorithm.

The following is a concise description of the algorithms used by the three fitting functions (note that the current version of the R package fitHeavyTail does not allow yet a regularization term with a target):


Aas, K., and Haff, I. H. (2006). The generalized hyperbolic skew student’st-distribution. Journal of Financial Econometrics, 4(2), 275–309.
Lee, S., and McLachlan, G. J. (2014). Finite mixtures of multivariate skew t-distributions: Some recent and new results. Statistics and Computing, 24(2), 181–202.
Liu, C., and Rubin, D. B. (1995). ML estimation of the t-distribution using EM and its extensions, ECM and ECME. Statistica Sinica, 5(1), 19–39.
Liu, C., Rubin, D. B., and Ying Nian Wu, and. (1998). Parameter expansion to accelerate EM: The PX-EM algorithm. Biometrika, 85(4), 755–770.
Ollila, E., Palomar, D. P., and Pascal, F. (2021). Shrinking the eigenvalues of M-estimators of covariance matrix. IEEE Transactions on Signal Processing, 69, 256–269.
Pascal, F., Ollila, E., and Palomar, D. P. (2021). Improved estimation of the degree of freedom parameter of multivariate t-distribution.
Sun, Y., Babu, P., and Palomar, D. P. (2014). Regularized Tyler’s scatter estimator: Existence, uniqueness, and algorithms. IEEE Trans. Signal Processing, 62(19), 5143–5156.
Sun, Y., Babu, P., and Palomar, D. P. (2015). Regularized robust estimation of mean and covariance matrix under heavy-tailed distributions. IEEE Trans. Signal Processing, 63(12), 3096–3109.
Sun, Y., Babu, P., and Palomar, D. P. (2017). Majorization-minimization algorithms in signal processing, communications, and machine learning. IEEE Transactions on Signal Processing, 65(3), 794–816.
Zhou, R., Liu, J., Kumar, S., and Palomar, D. P. (2019). Robust factor analysis parameter estimation. Lecture Notes in Computer Science (LNCS).