Background

The package MultiStatM provides general formulae for set partitions, multivariate moments and cumulants, vector Hermite polynomials. It provides theoretical formulae for some important symmetric and asymmetric multivariate distributions and well as estimation functions for multivariate moments and cumulants and connected measures of multivariate skewness and kurtosis.

The formulae implemented in the package can be found in the book “Multivariate Statistical Methods - Going Beyond the Linear”, Springer 2021 by Gy.Terdik and are fully general. For example, in the conversion formulae from multivariate moment to multivariate cumulants, given any list of (numerical) multivariate moments up to order \(k\), the conversion formula provides all multivariate cumulants up to order \(k\); this differs to a large degree from the formulae provided in the package kStatistics (Di Nardo and Guarino, 2022) which calculates one by one (individually) the cumulants of order \(r\) which are the entries of our cumulant vectors.

The packages MaxSkew and MultiSkew (Franceschini and Loperfido (2017a,b)) for detecting, measuring and removing multivariate skewness, computes the third multivariate cumulant of either the raw, centered or standardized data; s the main measures of multivariate skewness, together with their bootstrap distributions and provides orthogonal data projections with maximal skewness.

The package matrixcalc (Novomestky (2021)) provides the Commutation matrix, Elimination matrix, Duplication matrix for Cartesian tensor products of two vectors, which are particular cases of those provided in the package MultiStatM.

The package sn ( Azzalini (2022)) discusses for the skew-normal and the skew-t distributions, statistical methods are provided for data fitting and model diagnostics, in the univariate and the multivariate case. Random numbers generator for multivariate skew distributions are provided. In the package MultiStatM complete formulae for theoretical multivariate moments and cumulants of any order are implemented.

The package moments (Komsta and Novomestky (2022)) deals with functions to calculate moments, cumulants, Pearson’s kurtosis, Geary’s kurtosis and skewness; tests related to them from univariate data.

A careful study of the cumulants is a necessary and typical part of nonlinear statistics. Such a study of cumulants for multivariate distributions is made complicated by the index notations. One solution to this problem is the usage of tensor analysis. In this package (and the connected book) we offer an alternate method, which we believe is simpler to follow. The higher-order cumulants with the same degree for a multivariate vector can be collected together and kept as a vector. To be able to do so, we introduce a particular differential operator on a multivariate function, called the T -derivative, and use it to obtain cumulants and provide results which are somewhat analogous to well-known results in the univariate case.

More specifically, with the symbol \(\otimes\) denoting the Cartesian tensor product, consider the operator \(D_{\boldsymbol{\lambda}}^{\otimes}\), which we refer to as the \(\operatorname{T}\)-derivative; see Jammalamadaka et al. (2006) for details. For any function \(\boldsymbol{\phi}(\boldsymbol{\lambda})\), the~\(\operatorname{T}\)-derivative is defined as \[\begin{equation}\label{Tderiv} D_{\boldsymbol{\lambda}}^{\otimes}\boldsymbol{\boldsymbol{\phi}}% (\boldsymbol{\lambda})=\operatorname{vec}\left(\left( \frac{\partial\boldsymbol{\phi }(\boldsymbol{\lambda})}{\partial\boldsymbol{\lambda}^{\top}}\right) ^{\top }\right)=\boldsymbol{\phi}(\boldsymbol{\lambda})\otimes\frac{\partial}{\partial \boldsymbol{\lambda}}.% \end{equation}\] \({\boldsymbol{\phi}}\) is \(k\)-times differentiable, with~its \(k\)-th \(\operatorname{T}\)-derivative \(D_{\boldsymbol{\lambda}}^{\otimes k}\boldsymbol{\boldsymbol{\phi} }(\boldsymbol{\lambda})=D_{\boldsymbol{\lambda}}^{\otimes}\left( D_{\boldsymbol{\lambda}}^{\otimes k-1}\boldsymbol{\boldsymbol{\phi} }(\boldsymbol{\lambda})\right)\).

In the following we demonstrate the use of this technique through the characterization of several multivariate distributions via their cumulants and by extending the discussion to statistical inference for multivariate skewness and kurtosis.

We note that Kollo (2006) provides formulae for cumulants in terms of matrices; however, retaining a matrix structure for all higher-order cumulants leads to high-dimensional matrices with special symmetric structures which are quite hard to follow notionally and computationally. McCullagh (2018) provides quite an elegant approach using tensor methods; however, tensor methods are not very well known and computationally not so simple.

The method discussed here is based on relatively simple calculus. Although the tensor product of Euclidean vectors is not commutative, it has the advantage of permutation equivalence and allows one to obtain general results for cumulants and moments of any order, as it will be demonstrated in this paper, where general formulae, suitable for algorithmic implementation through a computer software, will be provided.

Methods based on a matrix approach do not provide this type of result; see also (Ould-Baba (2015), which goes as far as the sixth-order moment matrices, whereas there is no such limitation in our derivations and our results. For further discussion, one can see also Kolda (2009) and Qi (2006).

In MultiStatM 2.0.0 many functions of the previous version 1.2.1 have been renamed or grouped for better clarity and joining similar functions producing, for example the same output for univariate or multivariate cases

The table below provides a complete plan of transition from MultiStatM 1.2.1 to MultiStatM 2.0.0. Functions of version 1.2.1 which are within the same rowhave been grouped into a single function. For example the functions conv_Cum2Mom and conv_Cum2MomMulti which provided cumulants to moments conversion respectively in the univariate and multivariate cases have been joined in the function Cum2Mom which has now an option Type=c("Univariate","Multivariate").

Set Partitions

MultiStatM provides several functions dealing with set partitions. Such functions provide some basic tools used to build the multivariate formulae for moments and cumulants in the following sections.

Generally a set of \(N\) elements can be split into a set of disjoint subsets, i.e. it can be partitioned. The set of \(N\) elements will correspond to set \(1 : N = \{1, 2, \dots ,N\}\). If \({\cal{K}} = \{b_1, b_2, \dots , b_r \}\) where each \(b_j \subset 1 : N\), then \({\cal{K}}\) is a partition provided \(\cup b_j = 1 : N\), each \(b_j\) is non-empty and \(b_j \cap b_i = \emptyset\) (the empty set) is disjoint whenever \(j \neq i\). The subsets \(b_j\), \(j = 1, 2, \dots, r\) are called the blocks of \(\cal{K}\). We will call \(r\) (the number of the blocks in partition \(\cal{K}\)), the size of \(\cal{K}\), and denote it by \(|{\cal{K}}| = r\), and a partition with size \(r\) will be denoted by \({\cal{K}}_{\{r\}}\). Let us denote the set of all partitions of the numbers \(1 : N\) by \({\cal{P}}_N\).

Consider next a partition \({\mathcal{K}}_{\{r\}}=\{b_{1},b_{2},\dots,b_{r}\}\in {\mathcal{P}}_{N}\), with size \(r\). Denote the cardinality \(k_{j}\) of a block in the partition \({\mathcal{K}}_{\{r\}}\), i.e. \(k_{j}=|b_{j}|\). The type of a partition \({\mathcal{K}}_{\{r\}}\) is \(l=[l_{1},\dots ,l_{N}]\), if \({\mathcal{K}}_{\{r\}}\) contains exactly \(l_{j}\) blocks with cardinality \(j\). The type \(l\) is with length \(N\) always. A partition with size \(r\) and type \(l\) will be denoted by \({\mathcal{K}}_{\{r|l\}}\). It is clear that \(l_{j}\geq 0\), and \(\sum_{j}jl_{j}=N\), and \(\sum_{j}l_{j}=r\). Naturally, some \(l_{j}\)’s are zero. A block constitutes a row vector of entries \(0\)’s and \(1\)’s with length \(N\). The places of \(1\)’s correspond to the elements of the block. A partition matrix collects the rows of its blocks, it is an \(r\times N\) matrix with column-sums \(1\).

The basic function is PartitionTypeAll which provides complete information on the partition of a set of N elements, namely:

• S_N_r: a vector with the number of partitions of size r=1, r=2, etc. (Stirling numbers of the second kind); S_N_r[r] denotes the number of partition matrices of size r.

• Part.class: the list of all possible partitions given as partition matrices. This list is enumerated according to S_N_r[r], \(r=1,2,\ldots N\), such that the partition matrices with size R are listed from \(\sum_{r<R}\)\(\left[ r\right] +1\) up to \(\sum_{r\leq R}\)\(\left[ r\right]\). The order of the partition matrices within a fixed size is called canonical.

• S_r_j: a list of vectors of number of partitions with given types grouped by partitions of size r=1, r=2, etc.; an entry is the number of partitions with that type.

• eL_r: a list of partition types with respect to partitions of size r=1, r=2, etc. ; since a partition type is a row vector with length N this list includes matrices of types (row vectors), the number of rows are the length of vectors of S_r_j of a given size r.

Example 1. Consider the case where N=4 and run the following

PTA<-PartitionTypeAll(4)

S_N_r provides the number of partitions with r=1 to r=4 blocks:

PTA$S_N_r
#> [1] 1 7 6 1

All the partition matrices are listed in Part.class, for example the first partition matrix of size 2 among 7 is

PTA$Part.class[[2]]
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    1    1    0
#> [2,]    0    0    0    1

The second partition matrix of size 3 is

PTA$Part.class[[10]] 
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    0    1    0
#> [2,]    0    1    0    0
#> [3,]    0    0    0    1

etc.. If one interested in the number of partitions with different types for r=2 then consider the list S_r_j, i.e.

PTA$S_r_j[[2]] 
#> [1] 4 3

That is, for partitions with r=2 blocks, there are 2 possible types with 4 and 3 partitions each. These types will show up in the list eL_r:

PTA$eL_r
#> [[1]]
#> [1] 0 0 0 1
#> 
#> [[2]]
#>      [,1] [,2] [,3] [,4]
#> [1,]    1    0    1    0
#> [2,]    0    2    0    0
#> 
#> [[3]]
#> [1] 2 1 0 0
#> 
#> [[4]]
#> [1] 4 0 0 0

From the results above we see that there are: 1 partition of 1 block ( r=1); 7 partitions of two blocks (r=2); 6 partitions of 3 blocks and 1 partition of 4 blocks. From PTA$eL_r}[[2]] we see that there are two types of partition with r=2: the first is of type \(\left[ 1,0,1,0\right] ={(l_{1}=1,l_{2}=0,l_{3}=1,l_{4}=0)}\) with 4 partitions and the second is of type \(\left[ 0,2,0,0\right] ={(l_{1}=0,l_{2}=2,l_{3}=0,l_{4}=0)}\) with 3 partitions.

Another general function in this class is Partitions which has a Type argument in order to specify the type of partition to compute: 2Perm which provides the permutation of N elements according to a partition matrix L; Diagram which provides the list of partition matrices indecomposable with respect to L, representing diagrams without loops; Indecomp, which provides the list of all indecomposable partitions with respect to a partition matrix L; Pairs, which provides the list of partitions dividing into pairs the set of N elements.

Commutators, symmetrizer and selection matrices

The CommutatorMatr function produces commutators and selection matrices. The use of matrices allows represent as linear combinations problems of permutation and powers of T-products. On the other side, the size of these matrix can quickly become quite important. To deal with this issues and option for sparse matrices is always provided; also a corresponding CommutatorIndx function is provided; these function provide selection vectors which give equivalent results and the corresponding functions in the group Matr.

Kmn produces a commutation matrix, with usual notation \(\mathbf{K}_{m \cdot n}\), of dimension \(mn \times mn\) such that, given a matrix \(\mathbf{A}\) \(m\times n\), \(\mathbf{K}_{m \cdot n} \operatorname{vec}\mathbf{A}=\operatorname{vec}\mathbf{A}'\) (see @terdik2021multivariate, p.8) while Kperm produce any permutation of Kronecker products of vectors of any length.

Example 2. For the product of vectors \(\mathbf{a}_1 \otimes \mathbf{a}_2 \otimes\mathbf{a}_3\) of dimensions \(d_1\) to \(d_3\) respectively. CommutatorMatr(Type="Kperm",c(3,1,2),c(d1,d2,d3)) produces \(\mathbf{a}_3 \otimes \mathbf{a}_1 \otimes\mathbf{a}_2\).

a1<-c(1,2)
a2<-c(2,3,4)
a3<-c(1,3)
p1<-a1%x%a2%x%a3
c(CommutatorMatr(Type="Kperm",c(3,1,2),c(2,3,2))%*%p1) 
#>  [1]  2  3  4  4  6  8  6  9 12 12 18 24
a3%x%a1%x%a2 
#>  [1]  2  3  4  4  6  8  6  9 12 12 18 24

The same result can be obtained by using CommutatorIndx

p1[CommutatorIndx(Type="Kperm",c(3,1,2),c(2,3,2))]
#>  [1]  2  3  4  4  6  8  6  9 12 12 18 24

The CommutatorMatr with Type="Mixing" is exploited for deriving the covariance matrix of Hermite polynomials; see Terdik (2021, 4.6).

The Elimination and Qplication matrices- related functions respectively eliminate and restore duplicated or q-plicated elements in powers of T-products.

a<-c(1,2)
a3<-a%x%a%x%a
a3
#> [1] 1 2 2 4 2 4 4 8
c(EliminMatr(2,3)%*%a3)
#> [1] 1 2 4 8
c(QplicMatr(2,3)%*%EliminMatr(2,3)%*%a3)
#> [1] 1 2 2 4 2 4 4 8

Closely connected to the above matrices are the functions UnivMomCum and EliminIndx. The former provides a vector of indexes to select univariate moments or cumulants of the single elements of a d-vector X from available vector of T-moments and T-cumulants. The latter eliminates the duplicated/q-plicated elements in a T-vector of multivariate moments and cumulants. The function EliminIndx produces the same results as EliminMatr and it is less demanding in terms of memory. The use of EliminMatr can be preferable is one wishes to deal with linear combination of matrices. See examples 4 and 6 below for the use of UnivMomCum and EliminIndx.

The symmetrizer matrix, a \(d^n \times d^n\) matrix for the symmetrization of a T-product of \(n\) vectors with the same dimension \(d\) which overcomes the difficulties arising from the non commutative property of the Kronecker product, and simplifies considerably the computation formulae for multivariate polynomials and their derivatives (see Holmquist (1996) for details). The symmetrizer for a T-product of \(q\) vectors of dimension \(d\) is defined as \[ \mathbf{S}_{d \mathbf{1}q}=\frac{1}{q} \sum_{p \in \cal{P}_q} \mathbf{K}_p \] where \(\cal{P}_q\) is the set of all permutations of the numbers \(1:q\) and \(\mathbf{K}_p\) is the commutator matrix of for the permutation \(p \in \cal{P}_q\), (i.e. the CommutatorMatrKperm with Type="Kperm"of the package). Note that, by definition, computing the symmetrizer requires \(q!\) operations; in the package, the computational complexity is overcome by exploiting the Chacon and Duong (2015) efficient recursive algorithms for functionals based on higher order derivatives.

Multivariate T-Hermite Polynomials

Consider a Gaussian vector \(\mathbf{X}\) of dimension \(d\) with \(\operatorname{E}\mathbf{X}\) and \(\mathbf{\Sigma}=\operatorname{Cov}(\mathbf{X})=\operatorname{E}\mathbf{X X}'\) and define the generator function \[\begin{split} \Psi(\mathbf{X}; \mathbf{a})&=\exp \left(\mathbf{a}'\mathbf{X} - \frac{1}{2} \mathbf{a}' \mathbf{\Sigma} \mathbf{a}\right) \\ &=\exp \left(\mathbf{a}'\mathbf{X} - \frac{1}{2} \boldsymbol{\kappa}_2^{\otimes\prime} \mathbf{a}^{\otimes 2} \right) \\ \end{split} \] where \(\mathbf{a}\) is a \(d\)-vector of constants and \(\boldsymbol{\kappa}_2^{\otimes}=\operatorname{vec}\mathbf{\Sigma}\). The vector Hermite polynomials is defined via the T-derivative of the generator function, viz. \[ \mathbf{H}_n(\mathbf{X}) = D_\mathbf{a}^{\otimes n}\Psi(\mathbf{X};\mathbf{a})\big|_{\mathbf{a}=0} \] For example one has \[ \mathbf{H}_1(\mathbf{X})=\mathbf{X}, \quad \mathbf{H}_2(\mathbf{X})=\mathbf{X}^{\otimes 2} - \boldsymbol{\kappa}_2^{\otimes} \] Note that the multivariate T-Hermite polynomial \(\mathbf{H}_n(\mathbf{X})\) is a vector of dimension \(d^n\) which contains the \(n\)-th order polynomials of the vector \(\mathbf{X}^{\otimes n}\). For example the entries of \(\mathbf{H}_2(\mathbf{X})\) are the second order Hermite polynomials \(H_2(X_i,X_j)\), \(i,j=1,2, \dots d\); for \(d=2\) \[ \mathbf{H}_2(\mathbf{X}) = \left( (X_1^2 - \sigma_{11}), (X_1 X_2 - \sigma_{12}), (X_2 X_1 - \sigma_{21}), (X_2^2 - \sigma_{22})\right)^\prime. \] Note that \(\mathbf{H}_n(\mathbf{X})\) is \(n\)-symmetric, i.e. \(\mathbf{H}_2(\mathbf{X}) = \mathbf{S}_{d \mathbf{1}_n} \mathbf{H}_2(\mathbf{X})\) where \(\mathbf{S}_{d \mathbf{1}_n}\) is the symmetrizer defined in […]. From this one can get useful recursion formulae \[ \mathbf{H}_n(\mathbf{X})=\mathbf{S}_{d \mathbf{1}_n}\left( \mathbf{H}_{n-1}(\mathbf{X}) \otimes \mathbf{X}- (n-1) \mathbf{H}_{n-2}(\mathbf{X}) \otimes \boldsymbol{\kappa}_2^{\otimes} \right). \]

For further details, consult Terdik (2021, 4.5).

The definition of the \(d\)-variate Hermite polynomial requires the covariance matrix \(\mathbf{\Sigma}\) of the vector \(\mathbf{X}\). The HermiteN(...,Type="Univariate") and HermiteN(...,Type="Multivariate") functions compute the univariate and d-variate Hermite polynomials and their inverses up to a given order N evaluated at \(x\) for a given covariance matrix Sig2. By default Sig2=\(I_\mathbf{d}\).

Example 3 The first and the second \(3\)-variate Hermite polynomials evaluated at x<-c(1,2,3) where \(x\) is the realization of \(\mathbf{X} \sim N(\mathbf{0}, I_{\mathbf{3}})\) is

x<-c(1,2,3)
H2<-HermiteN(x,N=2,Type="Multivariate")
H2[[1]]
#> [1] 1 2 3
H2[[2]]
#> [1] 0 2 3 2 3 6 3 6 8

If x is the realization of \(\mathbf{X} \sim N(\mathbf{0}, 4I_\mathbf{2})\)

H2<-HermiteN(x,Sig2=4*diag(3),N=2,Type="Multivariate")
H2[[1]]
#> [1] 1 2 3
H2[[2]]
#> [1] -3  2  3  2  0  6  3  6  5

One can recover the vector x from H2 with the inverse function:

HermiteN2X(H2,N=2,Sig2=4*diag(3),Type="Multivariate")[[1]] 
#> [1] 1 2 3

The function HermiteCov12 can be exploited to obtain the covariance matrix of \(H_N(\mathbf{X}_1)\) and \(H_N(\mathbf{X}_2)\) for vectors \(\mathbf{X}_1\) and \(\mathbf{X}_2\) having covariance matrix \(\mathbf{\Sigma_{12}}\).

Covmat<-matrix(c(1,0.8,0.3,0.8,2,1,0.3,1,2),3,3)
Cov_X1_X2 <- HermiteCov12(SigX12=Covmat,N=3)

Multivariate T-moments and T-cumulants

Multivariate moments and cumulants of all orders of a random vector \(\mathbf{X}\) in \(d\)-dimensions, with mean vector \(\boldsymbol{\mu}\) and covariance matrix \(\mathbf{\Sigma}\) can be obtained by applying the T-derivative respectively to the characteristic function and the log of the CF.

More formally, let \(\boldsymbol{\lambda}\) a \(d\)-vector of real constants; \(\phi_{\mathbf{X}}(\boldsymbol{\lambda})\) and \(\psi_{\mathbf{X}}(\boldsymbol{\lambda})=\log\phi_{\mathbf{X}}(\boldsymbol{\lambda})\) denote, respectively, the characteristic function and the cumulant- function of \(\mathbf{X}\).

Then the \(k\)-th order moments and cumulants of the vector \(\mathbf{X}\) are obtained as \[ \boldsymbol{\mu}^\otimes_{\mathbf{X},k} = (-\mathbf{i})^k D_{\boldsymbol{\lambda}}^{\otimes k}\boldsymbol{\boldsymbol{\phi}% }_{\mathbf{X}}(\boldsymbol{\lambda}) \big|_{\boldsymbol{\lambda}=0}. \] \[ \boldsymbol{\kappa}^\otimes_{\mathbf{X},k} = \underline{\operatorname{Cum}}_k(\mathbf{X})= (-\mathbf{i})^k D_{\boldsymbol{\lambda}}^{\otimes k}\boldsymbol{\boldsymbol{\psi}% }_{\mathbf{X}}(\boldsymbol{\lambda}) \big|_{\boldsymbol{\lambda}=0}. \] Note that \(\boldsymbol{\mu}_{\mathbf{X},k} = \operatorname{E}\mathbf{X}^{\otimes k}\) that is a vector of dimension \(d^k\) that contains all possible moments of order order \(k\) formed by \(X_1, \dots, X_d\). This approach has the advantage of being straightforwardly extendable to any \(k\)-th order moment. An analogous discussion can be done for cumulants.

Note that one has \(\boldsymbol{\kappa}_{\mathbf{X},2} =\operatorname{vec} \mathbf{\Sigma}\).

The package MultiStatM contains functions which obtains moments from cumulants and vice-versa as well as function which provide theoretical moments and cumulants for some important multivariate distributions.

The Cum2Mom and Mom2Cum either for the univariate and multivariate cases provide conversion formulae for cumlants from moments and viceversa given any list of (theoretical) moments (or cumulants).

The conversion formula from moments to cumulants (see Terdik (2001, 3.4)) is given by \[\begin{split} \boldsymbol{\mu}_{n}^\otimes &= \sum_{\cal{K} \in \cal{P}_n} \mathbf{K}^{-1}_{p(\cal{K})} \prod^\otimes_{b_j \in \cal{K}} \kappa^\otimes_{|b_j|}\\ &= \mathbf{S}_{d \mathbf{1}_n}\left( \sum_{r=1}^n \sum_{\sum l_j =r, \sum j l_j = n} \frac{n!}{\prod_{j=1}^n l_j! (j!)^{l_j}} \prod_{j=1:n-r+1}^\otimes \kappa^{\otimes l_j}_j\right) \end{split} \] where the summation is over all partitions \(\cal{K} = \{b_1, b_2,\dots, b_k\}\) of \(1 : n\); \(|b_j|\) denotes the cardinality of block \(b_j\). The simpler second formula, exploiting the symmetrizer matrix, derives from symmetry of \(\boldsymbol{\mu}_{n}^\otimes\).

As far as the formula from cumulants to moments (Terdik (2021, 3.4)) is concerned, \[ \boldsymbol{\kappa}_{n}^\otimes = \mathbf{S}_{d \mathbf{1}_n}\left( \sum_{r=1}^n (-1)^{r-1} (r-1)!\sum_{\sum l_j =r, \sum j l_j = n} \prod_{j=1:n-r+1}^\otimes \frac{1}{{l_j}!}\left( \frac{1}{j!}\boldsymbol{\mu}^{\otimes}_j\right)^{l_j}\right) \]

Example 4. Consider the case of the 2-variate standard normal distribution with null mean vector and covariance matrix with unit elements on the main diagonal and off-diagonal elements equal to 0.5; in this case the the first four moments are given in the vector mu below

mu<-list(c(1,1),c(2,1.5,1.5,2),c(4,3,3,3,3,3,3,4),c(10,7,7,6.5,7,6.5,6.5,7,7,6.5,6.5,7,6.5,7,7,10))
cum<-Mom2Cum(mu, Type="Multivariate")
cum
#> [[1]]
#> [1] 1 1
#> 
#> [[2]]
#> [1] 1.0 0.5 0.5 1.0
#> 
#> [[3]]
#> [1] 0 0 0 0 0 0 0 0
#> 
#> [[4]]
#>  [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Getting back to moments

Cum2Mom(cum,Type="Multivariate")
#> [[1]]
#> [1] 1 1
#> 
#> [[2]]
#> [1] 2.0 1.5 1.5 2.0
#> 
#> [[3]]
#> [1] 4 3 3 3 3 3 3 4
#> 
#> [[4]]
#>  [1] 10.0  7.0  7.0  6.5  7.0  6.5  6.5  7.0  7.0  6.5  6.5  7.0  6.5  7.0  7.0
#> [16] 10.0

I one wishes to select only the distinct moments from the vector of third moments, then

mu[[3]][EliminIndx(2,3)]
#> [1] 4 3 3 4

Alternatively one can also use the elimination matrix

r.mu<-EliminMatr(2,3)%*% mu[[3]]
c(r.mu)
#> [1] 4 3 3 4

Note that EliminMatr does not correspond the the function unique, rather it individuates the duplicated elements from the symmetry of the Kronecker product. This allow to recover the whole vector when needed.

c(QplicMatr(2,3)%*%r.mu)
#> [1] 4 3 3 3 3 3 3 4

The same result by using QplicIndx

r.mu[QplicIndx(2,3)]
#> [1] 4 3 3 3 3 3 3 4

The MomCum functions provide theoretical moments and cumulants for some common multivariate distributions: Skew-normal, Canonical Fundamental Skew-normal (CFUSN), Uniform distribution on the sphere, central folded Normal distribution (univariate and multivariate); for detail on the multivariate formulae used see @jamma2021San. Evaluation of theoretical moments and cumulants is done by the MomCumNAME group of functions. Some more details on the multivariate distributions considered are reported in the list below.

• A \(d\)-vector \(\mathbf{U}\) having uniform distribution on the sphere \(\mathbb{S}_{d-1}\). Moments and cumulants of all orders are provided for \(\mathbf{U}\) by the function MomCumUniS; the function EVSKUniS can compute moments and cumulants (up to the 4th order), skewness, and kurtosis of \(\mathbf{U}\) (Type="Standard") and its modulus (Type="Modulus"). Recall that any \(d\)-vector, say \(\mathbf{W}\), has a spherically symmetric distribution if that distribution is invariant under the group of rotations in \(\mathbb{R}^{d}\). This is equivalent to saying that \(\mathbf{W}\) has the stochastic representation \(\mathbf{W}=R\mathbf{U}\) where \(R\) is a non negative random variable. Moments and cumulants of \(\mathbf{W}\) can be obtained by its stochastic representation as discussed in @jamma2021San, Theorem 1 and Lemma 1. Furthermore a \(d\)-vector \(\mathbf{X}\) has an elliptically symmetric distribution if it has the representation \[ \mathbf{X}=\boldsymbol{\mu}+\boldsymbol{\Sigma}^{1/2}\mathbf{W}% \] where \(\boldsymbol{\mu}\in\mathbb{R}^{d}\), \(\boldsymbol{\Sigma}\) is a variance-covariance matrix and \(\mathbf{W}\) is spherically distributed. Hence the cumulants of \(\mathbf{X}\) are just constant times the cumulants of \(\mathbf{W}\) except for the mean i.e. \[ \underline{\operatorname*{Cum}}_{m}\left( \mathbf{X}\right) =\left( \boldsymbol{\Sigma}^{1/2}\right) ^{\otimes m}\underline{\operatorname*{Cum}% }_{m}\left( \mathbf{W}\right) . \]

• If \(\mathbf{Z}\) denotes a \(d\)-vector with \(d\)-variate standard normal distribution, the function MomCumZabs provide the moments and cumulants of \(|\mathbf{Z}|\) respectively in the univariate (Type="univariate") and multivariate case (Type="Multivariate").

• The multivariate skew-normal distribution introduced by @azzalini1996multivariate, whose marginal densities are scalar skew-normals. A \(d\)-dimensional random vector \(\mathbf{X}\) is said to have a multivariate skew-normal distribution, \(\text{SN}_{d}\left(\boldsymbol{\mu},\boldsymbol{\Omega},\boldsymbol{\alpha}\right)\) with shape parameter \(\boldsymbol{\alpha}\) if it has the density function \[ 2\varphi\left( \mathbf{X};\boldsymbol{\mu},\boldsymbol{\Omega}\right) \Phi\left( \boldsymbol{\alpha}^{\top}\left( \mathbf{X}-\boldsymbol{\mu }\right) \right) , \quad\mathbf{X} \in\mathbb{R}^{d}, \] where \(\varphi\left(\mathbf{X};\boldsymbol{\mu},\boldsymbol{\Omega}\right)\) is the \(d\)-dimensional normal density with mean \(\boldsymbol{\mu}\) and correlation matrix \(\boldsymbol{\Omega}\); here \(\varphi\) and \(\Phi\) denote the univariate standard normal density and the cdf. For a general formula for cumulants, see @jamma2021San, Lemma 4. For this distribution are available the functions MomCumSkewNorm, which computes the theoretical values of moments and cumulants up to the r-th order and EVSKSkewNorm which gives mean vector, covariance, skewness and kurtosis vectors and other measures.

• @arellano2005fundamental introduced the CFUSN distribution (cf. their Proposition 2.3), to include all existing definitions of skew-normal distributions. The marginal stochastic representation of \(\mathbf{X}\) with distribution \(\text{CFUSN}_{d,m}\left(\boldsymbol{\Delta}\right)\) is given by \[ \mathbf{X}=\boldsymbol{\Delta}\left\vert \mathbf{Z}_{1}\right\vert +\left( \mathbf{I}_{d}-\boldsymbol{\Delta\Delta}^{\top}\right) ^{1/2}\mathbf{Z}_{2} \] where \(\boldsymbol{\Delta}\), is the \(d\times m\) skewness matrix such that \(\left\Vert \boldsymbol{\Delta}\underline{a}\right\Vert <1\), for all \(\left\Vert \underline{a}\right\Vert =1\), and \(\mathbf{Z}_{1}\in\mathcal{N}\left( 0,\mathbf{I}_{m}\right)\) and \(\mathbf{Z}_{2}\in\mathcal{N}\left( 0,\mathbf{I}_{d}\right)\) are independent (Proposition 2.2. Arellano-Valle and Genton (2005)). A simple construction of \(\boldsymbol{\Delta}\) is \(\boldsymbol{\Delta}=\boldsymbol{\Lambda}\left(\mathbf{I}_{m}\mathbf{+}\boldsymbol{\Lambda}^{\top}\boldsymbol{\Lambda}\right)^{-1/2}\) with some real matrix \(\boldsymbol{\Lambda}\) with dimensions \(d\times m\). The \(\text{CFUSN}_{d,m}\left(\boldsymbol{\mu},\boldsymbol{\Sigma},\boldsymbol{\Delta}\right)\) can be defined via the linear transformation \(\boldsymbol{\mu}+\boldsymbol{\Sigma}^{1/2}\mathbf{X}\). For a general formula for cumulants, see @jamma2021San, Lemma 5. For this distributions MomCumCFUSN provides moments and cumulants up to the \(r\)-th order.

The Random generation family of functions in provide random number generators for multivariate distributions.

Example 5. For a skew-normal distribution with \(\alpha=(10,5,0)\) and correlation function \(\Omega= \text{diag} (1,1,1)\) we have the third moments and cumulants are

alpha<-c(10,5,0)
omega<-diag(3)
MSN<-MomCumSkewNorm(r=3,omega,alpha,nMu=TRUE)
round(MSN$Mu[[3]],3)
#>  [1] 1.568 0.073 0.000 0.073 0.570 0.000 0.000 0.000 0.711 0.073 0.570 0.000
#> [13] 0.570 0.996 0.000 0.000 0.000 0.355 0.000 0.000 0.711 0.000 0.000 0.355
#> [25] 0.711 0.355 0.000
round(MSN$CumX[[3]],3)
#>  [1] 0.154 0.077 0.000 0.077 0.039 0.000 0.000 0.000 0.000 0.077 0.039 0.000
#> [13] 0.039 0.019 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000
#> [25] 0.000 0.000 0.000

As another example, for the Uniform distribution on the sphere, the fourth cumulant is:

EVSKUniS(3,  Type="Standard")$Kurt.U
#>  [1] -1.2  0.0  0.0  0.0 -0.4  0.0  0.0  0.0 -0.4  0.0 -0.4  0.0 -0.4  0.0  0.0
#> [16]  0.0  0.0  0.0  0.0  0.0 -0.4  0.0  0.0  0.0 -0.4  0.0  0.0  0.0 -0.4  0.0
#> [31] -0.4  0.0  0.0  0.0  0.0  0.0 -0.4  0.0  0.0  0.0 -1.2  0.0  0.0  0.0 -0.4
#> [46]  0.0  0.0  0.0  0.0  0.0 -0.4  0.0 -0.4  0.0  0.0  0.0 -0.4  0.0  0.0  0.0
#> [61] -0.4  0.0  0.0  0.0  0.0  0.0  0.0  0.0 -0.4  0.0 -0.4  0.0 -0.4  0.0  0.0
#> [76]  0.0 -0.4  0.0  0.0  0.0 -1.2

Estimation

Estimating functions starting from a vector of multivariate data are available: multivariate moments and cumulants, skewness and kurtosis vectors Mardia’s skewness and kurtosis indexes, Mori, Rohatgi, Szekely (MRSz’s) skewness vector and kurtosis matrices.

A complete picture of skewness is provided by the third-order T-cumulant (skewness vector) of a standardized \(\mathbf{X}\); set \(\mathbf{Y}=\mathbf{\Sigma}^{-1/2}(\mathbf{X}-\boldsymbol{\mu})\), then the skewness vector is \[ \boldsymbol{\kappa}_{\mathbf{Y},3}^\otimes =\underline{\operatorname{Cum}}_3 \left( \mathbf{Y}\right)=\left(\mathbf{\Sigma}^{-1/2}\right)^{\otimes 3} \boldsymbol{\kappa}_{\mathbf{X},3}^\otimes. \] The total skewness of \(\mathbf{X}\) is defined by the square norm of the skewness vector: \(\gamma_{1,d}=\|\boldsymbol{\kappa}_{\mathbf{Y},3}^\otimes\|^2\). This definition guarantees that skewness is invariant under the shifting and orthogonal transformations, in other words it is affine invariant.

We note that Mardia’s multivariate skewness index (Mardia (1970)), denote it by \(\beta_{1,d}\), coincides with the total skewness \(\gamma_{1,d}\) since the third-order central moments and third-order cumulants are equal.

Mori, Rohatgi, Szekely (MRSz’s) skewness vector (Mori et al. (1994)) can also be recovered from the skewness vector as \[ \mathbf{b}(\mathbf{Y})= \left( \operatorname{vec}' \mathbf{I}_d \otimes \mathbf{I}_d \right)\boldsymbol{\kappa}_{\mathbf{Y},3}^\otimes \] Note that \(\operatorname{vec}' \mathbf{I}_d \otimes \mathbf{I}_d\) is a matrix of dimension \(d \times d^3\), which contains \(d\) unit values per-row, whereas all the others are 0; as a consequence, this measure does not take into account the contribution of cumulants of the type \(\operatorname{Cum}_3 (X_j,X_k,X_l)\), where all the three indices \(j\), \(k\), \(l\) are different from each other. The corresponding scalar measure of multivariate skewness is \(b(\mathbf{Y}) = \| \mathbf{b}(\mathbf{Y}) \|^2\).

The fourth-order T-cumulant of the standardized \(\mathbf{X}\), i.e. \(\boldsymbol{\kappa}_{\mathbf{Y},4}^\otimes\), will be called kurtosis vector of \(\mathbf{X}\); its square norm will be called the total kurtosis of \(\mathbf{X}\) \[ \gamma_{2,d}=\| \boldsymbol{\kappa}_{\mathbf{Y},4}^\otimes \|^2 \] Mardia’s kurtosis index \(\beta_{2,d}= \operatorname{E}\left( \mathbf{Y}'\mathbf{Y} \right)^2\) is related to the kurtosis vector by the formula \[ \beta_{2,d}= \left( \operatorname{vec}' \mathbf{I}_{d^2} \right)\boldsymbol{\kappa}_{\mathbf{Y},4}^\otimes +d(d+2) \] A consequence of this is that Mardia’s measure does not depend on all the entries of \(\boldsymbol{\kappa}_{\mathbf{Y},4}^\otimes\) which has \(d(d +1)(d +2)(d +3)/24\) distinct elements, while \(\beta_{2,d}\) includes only \(d^2\) elements among them. We note that if \(\mathbf{X}\) is Gaussian, then \(\boldsymbol{\kappa}_{\mathbf{Y},4}^\otimes=\mathbf{0}\).

Cardoso, Mori, Szekely, Rothagi define what we will call the CMRS kurtosis matrix \[ \mathbf{B}(Y) =\operatorname{E}\left( \mathbf{Y}\mathbf{Y}' \mathbf{Y}\mathbf{Y}' \right) -(d+2)\mathbf{I}_d \] which can be expressed in terms of the kurtosis vector as \[ \operatorname{vec}\mathbf{B}(Y)\left( \mathbf{I}_{d^2}\otimes \operatorname{vec}' \mathbf{I}_d \right)\boldsymbol{\kappa}_{\mathbf{Y},4}^\otimes \] Note also that \(\operatorname{tr} \mathbf{B}(Y) = \beta_{2,d}\).

For further discussion on the above indexes and further multivariate indexes of skewness and kurtosis, as well as their asymptotic theory one can consult Terdik (2021, section 6) and Jammalamadaka et al. (2021a,b).

The function Esti_Variance_Skew_Kurt provides estimates of the covariance matrix of the data-estimated skewness and kurtosis vectors (Terdik (2021), formulae 6.13 and 6.22).

Example 6. Consider a multivariate data vector of dimension \(d=3\) and \(n=250\) from the multivariate skew-normal distribution of Example 5. The estimated first four cumulants are listed in the object EsMSN obtained by the SampleEVSK function; the corresponding theoretical values are in the object ThMSN obtained by the EVSKSkewNorm function.

data<-rSkewNorm(1000,omega,alpha)
EsMSN<-SampleEVSK(data)
ThMSN<-EVSKSkewNorm(omega,alpha)

Compare the distinct elements of the estimated skewness vector and the theoretical ones using Eliminindx.

EsMSN$estSkew[EliminIndx(3,3)]   
#>  [1]  0.553645757  0.321903151 -0.042009525  0.214214390  0.007298792
#>  [6]  0.023751673  0.029318273  0.068464555 -0.049131406  0.065397721
ThMSN$SkewX[EliminIndx(3,3)]   
#>  [1] 0.68927167 0.34463583 0.00000000 0.17231792 0.00000000 0.00000000
#>  [7] 0.08615896 0.00000000 0.00000000 0.00000000

If one wishes to recover the estimated univariate skewness and kurtosis of the components \(X1\), \(X2\) and \(X3\) of \(X\), then, using UnivMomCum,

EsMSN$estSkew[UnivMomCum(3,3)]  ## Get univariate skewness for X1,X2,X3
#> [1] 0.55364576 0.02931827 0.06539772
EsMSN$estKurt[UnivMomCum(3,4)]  ## Get univariate kurtosis for X1,X2,X3
#> [1]  0.20291995 -0.03904306 -0.02658244

An estimate of Mardia’s skewness index is provided together with the p-value under the null hypothesis of normality. The theoretical value of Mardia’s skewness can be recovered from the element SkewX.tot in the object ThMSN.

SampleSkew(data,Type="Mardia")
#> $Mardia.Skewness
#> [1] 0.7887987
#> 
#> $p.value
#> [1] 2.345182e-23
ThMSN$SkewX.tot
#> [1] 0.9279208

The MRS skewness vector and index are provided together with the p-value for the skewness index under the null hypothesis of normality, The theoretical value, for the distribution at hand, can be computed using formula […]

SampleSkew(data,Type="MRSz")
#> $MRSz.Skewness.Vector
#> [1] 0.79161182 0.30209002 0.09185275
#> 
#> $MRSz.Skewness.Index
#> [1] 0.7263446
#> 
#> $p.value
#> [1] 1.164125e-15
as.vector(t(c(diag(3))%x%diag(3))%*%ThMSN$SkewX)
#> [1] 0.8615896 0.4307948 0.0000000

c(t(c(diag(3))%x%diag(3))%*%ThMSN$SkewX)  ## Theoretical MRS skewness vector
#> [1] 0.8615896 0.4307948 0.0000000

Acknowledgement

This work has been partially supported by the project TKP2021-NKTA of the University of Debrecen, Hungary.

References

Arellano-Valle, R.B., Genton, M.G. (2005) On fundamental skew distributions. Journal of Multivariate Analysis 96(1), 93–116.

Azzalini, A. (2022). The R package ‘sn’: The Skew-Normal and Related Distributions such as the Skew-t and the SUN (version 2.0.2). URL http://azzalini.stat.unipd.it/SN/, https://cran.r-project.org/package=sn

Azzalini, A,, Dalla Valle, A. (1996) The multivariate skew-normal distribution. Biometrika 83(4), 715–726

Chacón, J. E., & Duong, T. (2015). Efficient recursive algorithms for functionals based on higher order derivatives of the multivariate Gaussian density. Statistics and Computing, 25(5), 959-974.

Di Nardo, E. and Guarino, G. (2022). kStatistics: Unbiased Estimators for Cumulant Products and Faa Di Bruno’s Formula. R package version 2.1.1. https://CRAN.R-project.org/package=kStatistics

Franceschini, C. and Loperfido, N. (2017a). MultiSkew: Measures, Tests and Removes Multivariate Skewness. R package version 1.1.1. https://CRAN.R-project.org/package=MultiSkew

Franceschini, C. and Loperfido, N. (2017b). MaxSkew: Orthogonal Data Projections with Maximal Skewness. R package version 1.1. https://CRAN.R-project.org/package=MaxSkew

Holmquist, B. (1996). The d-variate vector Hermite polynomial of order. Linear Algebra and Its Applications, 237/238, 155–190.

Jammalamadaka, S. R., Subba Rao, T. and Terdik, Gy. (2006). Higher order cumulants of random vectors and applications to statistical inference and time series. Sankhya A, 68, 326–356.

Jammalamadaka, S. R., Taufer, E. & Terdik, Gy. H. (2021a). On multivariate skewness and kurtosis. Sankhya A, 83(2), 607-644.

Jammalamadaka, S. R., Taufer, E. & Terdik, Gy. H. (2021b). Asymptotic theory for statistics based on cumulant vectors with applications. Scandinavian Journal of Statistics, 48(2), 708-728.

Jammalamadaka,S.R. , Taufer,E. & Terdik, Gy. (2021c). Cumulants of Multivariate Symmetric and Skew Symmetric Distributions, Symmetry 13, 1383.

Kolda, T.G.; Bader, B.W. (2009). Tensor decompositions and applications. SIAM Rev. 51, 455–500.

Kollo, T. (2008). Multivariate skewness and kurtosis measures with an application in ICA. Journal of Multivariate Analysis 99(10), 2328–2338.

Komsta, L. and Novomestky F. (2022). moments: Moments, Cumulants, Skewness, Kurtosis and Related Tests. R package version 0.14.1. https://CRAN.R-project.org/package=moments

Novomestky, F. (2021). matrixcalc: Collection of Functions for Matrix Calculations. R package version 1.0-5. https://CRAN.R-project.org/package=matrixcalc

Qi, L. (2006). Rank and eigenvalues of a supersymmetric tensor, the multivariate homogeneous polynomial and the algebraic hypersurface it defines. J. Symb. Comput. 41, 1309–1327.

Mardia, K. V. (1970). Measures of multivariate skewness and kurtosis with applications. Biometrika, 57, 519–530.

McCullagh, P. (2018). Tensor methods in statistics. Chapman and Hall/CRC.

Móri, T. F., Rohatgi, V. K., & Székely, G. J. (1994). On multivariate skewness and kurtosis. Theory of Probability & Its Applications, 38(3), 547–551.

Ould-Baba, H.; Robin, V.; Antoni, J. (2015). Concise formulae for the cumulant matrices of a random vector. Linear Algebra Appl. 485, 392–416.

Terdik, Gy. (2021). Multivariate statistical methods - going beyond the linear. Springer.