The Mardia's multivariate skewness formula (Mardia, 1970) is
$$ b_{1, d} = \frac{1}{n^2}\sum^n_{i=1}\sum^n_{j=1}\left[
\left(\bold{X}_i - \bold{\bar{X}} \right)^{'} \bold{S}^{-1}
\left(\bold{X}_j - \bold{\bar{X}} \right) \right]^3, $$
where \(d\) is the number of variables, \(X\) is the target dataset
with multiple variables, \(n\) is the sample size, \(\bold{S}\) is
the sample covariance matrix of the target dataset, and \(\bold{\bar{X}}\)
is the mean vectors of the target dataset binded in \(n\) rows.
When the population multivariate skewness is normal, the
\(\frac{n}{6}b_{1,d}\) is asymptotically distributed as \(\chi^2\)
distribution with \(d(d + 1)(d + 2)/6\) degrees of freedom.