Given a \(p\)-dimensional multivariate random vector with mean vector \(\boldsymbol{\mu}\) and positive definite variance-covariance matrix \(\boldsymbol{\Sigma}\), Mardia's multivariate skewness is defined as
$$\beta_{1,p} = E[(\boldsymbol{X}_1 - \boldsymbol{\mu})' \boldsymbol{\Sigma}^{-1} (\boldsymbol{X}_2 - \boldsymbol{\mu})]^3,$$
where \(\boldsymbol{X}_1\) and \(\boldsymbol{X}_2\) are independently and identically distributed copies of \(\boldsymbol{X}\). For a multivariate random sample of size \(n\), \(\boldsymbol{x}_1, \boldsymbol{x}_1, \ldots, \boldsymbol{x}_n\), its sample version is defined as
$$\hat{\beta}_{1,p} = \frac{1}{n^2} \sum_{i=1}^{n} \sum_{j=1}^{n} [(\boldsymbol{x}_i - \bar{\boldsymbol{x}})'\boldsymbol{S}^{-1} (\boldsymbol{x}_j - \bar{\boldsymbol{x}})]^3,$$
where the sample mean \(\bar{\boldsymbol{x}} = \frac{1}{n}\sum_{i=1}^{n} \boldsymbol{x}_i\) and the sample variance-covariance matrix \(\boldsymbol{S} = \frac{1}{n} \sum_{i=1}^{n} (\boldsymbol{x}_i - \bar{\boldsymbol{x}}) (\boldsymbol{x}_i - \bar{\boldsymbol{x}})'\). It is assumed that \(n \ge p\).