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).
The package MultiStatM provides several functions
dealing with set partitions. Such functions provide some basic tools
used to built 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 \(\cal{K}_{\{r\}} = \{b_1, b_2, \dots , b_r \} \in \cal{P}_N\), with size \(r\). Denote the cardinality \(k_j\) of a block in the partition \(\cal{K}_{\{r\}}\), i.e. \(k_j =|b_j|\). The type of a partition \(\cal{K}_{\{r\}}\) is \(l = [ l_1, \dots , l_N]\), if \(\cal{K}_{\{r\}}\) contains exactly \(l_j\) blocks with cardinality \(j\). A partition with size \(r\) and type \(l\) will be denoted by \(\cal{K}_{\{r|l\}}\). It is clear that \(l_j ≥ 0\), and \(\sum_j jl_j=N\), and \(\sum_j l_j=r\). Naturally, some \(l_j\) ’s are zero.
| Function | Description |
|---|---|
| Partition_Type_All | Type and number of partitions |
| Partition_Indecomposable | Building indecomposable partitions |
| Partition_Pairs | Partition into pairs of the set 1:N |
| Permutation_Inverse | Inverse of a Permutation |
| Partition_DiagramsClosedNoLoops | Closed Diagrams without Loops |
The basic function is Partition_Type_All which provides
complete information on the partition of a set of N
elements, namely
Part.class: the list of all possible partitions
given as partition matrices
S_N_r: a vector with the number of partitions of
size r=1, r=2, etc. (Stirling numbers of second kind )
eL_r: a list of partition types with respect to
partitions of size r=1, r=2, etc.
S_r_j: vectors of number of partitions with given
types grouped by partitions of size r=1, r=2, etc.
Example 1. Consider the case where N=4
and run the following
All the partition matrices are listed in Part.class. See
the first three below
PTA$Part.class[[1]]
#> [1] 1 1 1 1
PTA$Part.class[[2]]
#> [,1] [,2] [,3] [,4]
#> [1,] 1 1 1 0
#> [2,] 0 0 0 1
PTA$Part.class[[3]]
#> [,1] [,2] [,3] [,4]
#> [1,] 1 1 0 1
#> [2,] 0 0 1 0The lists S_N_r and eL_r contain
respectively the number of partitions with r=1 blocks,
r=2 blocks, etc. and the types of partitions for each
partition of size r
PTA$S_N_r
#> [1] 1 7 6 1
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 0From 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 \((l_1=1,l_3=1)\) and the second is of type
\((l_2=2)\).
If one needs to know how many partitions of type \((l_1,l_3)\) are there for r=2
then consider the list S_r_j, i.e.
That is, for partitions with r=2 blocks, 4 are of type
\((l_1=1,l_3=1)\) and 3 are of type
\((l_2=2)\).
Other available partition functions are:
Partition_2perm which provides the permutation of N
elements according to a partition matrix \(pU\);
Partition_Indecomposable, which provides the list of
all indecomposable partitions with respect to a partition matrix
L
Partition_Pairs, which provides the list of
partitions dividing into pairs the set of N elements.
Partition_DiagramsClosedNoLoops Closed Diagrams
without Loops
Permutation_Inverse which provides the inverse of a
permutation
For further details see Terdik (2021, chapter 1).
| Function | Description |
|---|---|
| matr_Commutator_Kmn | Commutation matrix |
| matr_Commutator_Kperm | Commutator for T-products of vectors |
| matr_Commutator_Mixing | Mixing commutator |
| matr_Elimination | Elimination matrix |
| matr_Qplication | Qplication matrix |
| matr_Symmetry | Symmetrizer matrix |
| indx_Commutator_Kmn | Index vector for commutation of T-products of two vectors |
| indx_Commutator_Kperm | Index vector for commutation of T-products of any number of vectors |
| indx_Elimination | Distinct values selection vector |
| indx_Qplication | Qplication vector |
| indx_Symmetry | Symmetrizing vector |
| indx_UnivMomCum | Univariate moments and cumulants from T-vectors |
The matr group of functions produce 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 indx group of functions is
provided; these function provide selection vectors which give equivalent
results ans the corresponding functions in the group
matr.
The function matr_Commutator_Kmn produces a commutation
matrix, \(\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 Terdik (2021, p.8,(1.12)) ). The same result can be obtained using
indx_Commutator_Kmn
matr_Commutator_Kperm and
indx_Commutator_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.
matr_Commutator_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
as.vector(matr_Commutator_Kperm(c(3,1,2),c(2,3,2))%*%p1) ## same result as below
#> [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 24The same result can be obtained by using
indx_Commutator_Kperm
The matr_Commutator_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
as.vector(matr_Elimination(2,3)%*%a3)
#> [1] 1 2 4 8
as.vector(matr_Qplication(2,3)%*%matr_Elimination(2,3)%*%a3)
#> [1] 1 2 2 4 2 4 4 8Closely connected to the above matrices are the functions
indx_UnivMomCum and indx_Elimination. 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 indx_Elimination produces the same
results as matr_Elimination and it is less demanding in
terms of memory. The use of matr_Elimination can be
preferable is one wishes to deal with linear combination of matrices.
See examples 4 and 6 below for the use of indx_UnivMomCum
and indx_Elimination.
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 matr_Commutator_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.
| Function | Description |
|---|---|
| Hermite_Poly_HN | Univariate Hermite polynomials |
| Hermite_Poly_HN_Multi | Multivariate T-Hermite polynomials |
| Hermite_Poly_NH_Inv | Inverse univariate Hermite polynomial |
| Hermite_Poly_NH_Multi_Inv | Inverse of d-variate T-Hermite Polynomial |
| Hermite_N_Cov_X1_X2 | Computation of the covariance matrix between \(d\)-variate T-Hermite polynomials \(H_N(X_1)\) and \(H_N(X_2)\) |
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
Hermite_Poly_HN and Hermite_Poly_HN_multi
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
If x is the realization of \(\mathbf{X} \sim N(\mathbf{0},
4I_\mathbf{2})\)
H2<-Hermite_Poly_HN_Multi(x,Sig2=4*diag(3),N=2)
H2[[1]]
#> [1] 1 2 3
H2[[2]]
#> [1] -3 2 3 2 0 6 3 6 5One can recover the vector x from H2 with the inverse function:
The function Hermite_N_Cov_X1_X2 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}}\).
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 conv_Cum2Mom and conv_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).
| Function | Description |
|---|---|
| conv_Cum2Mom | Convert cumulants to moments (univariate) |
| conv_Cum2MomMulti | Convert T-cumulants to T-moments (multivariate) |
| conv_Mom2Cum | Convert moments to cumulants (univariate) |
| conv_Mom2CumMulti | Convert T-moments to T-cumulants (multivariate) |
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<-conv_Mom2CumMulti(mu)
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 0Getting back to moments
conv_Cum2MomMulti(cum)
#> [[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.0I one wishes to select only the distinct moments from the vector of third moments, then
Alternatively one can also use the elimination matrix
Note that matr_Elimination 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.
The same result by using indx_Qplication
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
(Jammalamadaka et al., 2021).
| Function | Description |
|---|---|
| distr_CFUSN_MomCum_Th | Moments and cumulants CFUSN |
| distr_ZabsM_MomCum_Th | Moments and cumulants multivariate central folded Normal distribution |
| distr_SkewNorm_MomCum_Th | Moments and cumulants d-variate Skew Normal |
| distr_Uni_MomCum_Th | Moments and cumulants Uniform Distribution on the Sphere |
| distr_Zabs_MomCum_Th | Moments and cumulants Central folded Normal distribution |
| distr_SkewNorm_EVSK_Th | EVSK multivariate Skew Normal |
| distr_UniAbs_EVSK_Th | Moments of the modulus of the Uniform distribution on the sphere |
| distr_Uni_EVSK_Th | EVSK Uniform on the sphere |
Expressions for theoretical moments and cumulants are provided by the
distr_NAME_MomCum_Th functions. In particular
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}\) (NAME=Uni) and
its modulus (NAME=UniAbs). 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
Jammalamadaka et al.(2021a, Theorem 1) and Jammalamadaka et al. (2021c,
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 Normal distribution, the
functions with NAME=Zabs and NAME=ZabsM
provide the moments and cumulants of \(|\mathbf{Z}|\) respectively in the
univariate (\(d=1\)) and multivariate
case.
The Multivariate Skew Normal distribution introduced by Azzalini and Dalla Valle (1996), 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 Jammalamadaka et al. (2021a, Lemma 4).
Arellano-Valle and Genton (2005) introduced the CFUSN distribution (cf. their Proposition 2.3), to include all existing definitions of SN 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 \(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 Jammalamadaka et al. (2021a, Lemma 5).
The Rand functions provide random number generators for
multivariate distributions.
The EVSK functions compute the theoretical values of the
mean vector, covariance, skewness vector, total skenwness, kurtosis
vector and total kurtosis for given multivariate distributions: Uniform
on the sphere, modulus of the Uniform distribution on the sphere,
Skew-normal distribution.
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<-distr_SkewNorm_MomCum_Th(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.000As another example, for the modulus of the Uniform distribution on the sphere, the fourth cumulant is:
distr_Uni_EVSK_Th(3, nCum = TRUE)$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.2Estimating 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.
| Function | Description |
|---|---|
| Esti_Hermite_Poly_HN_Multi | Estimation of the N-th d-variate polynomial |
| Esti_EVSK | Estimation of multivariate Mean, Variance, T-Skewness and T-Kurtosis vectors |
| Esti_Kurt_CMRSz | Estimation of Cardoso, Mori,Rohatgi , Szekely (CMRSz’s) kurtosis matrix |
| Esti_Kurt_Mardia | Estimation of Mardia’s Kurtosis Index |
| Esti_Kurt_Total | Estimation of the Total Kurtosis Index |
| Esti_MMom_MCum | Estimation of multivariate T-Moments and T-Cumulants |
| Esti_Skew_Mardia | Estimation of Mardia’s Skewness index |
| Esti_Skew_MRSz | Estimation for Mori, Rohatgi, Szekely (MRSz’s) skewness vector |
| Esti_Variance_Skew_Kurt | Estimated Variance of skewness and kurtosis vectors |
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 Esti_EVSK
function; the corresponding theoretical values are in the object
ThMSN obtained by the istr_SkewNorm_EVSK_Th
function.
data<-distr_SkewNorm_Rand(1000,omega,alpha)
EsMSN<-Esti_EVSK(data)
ThMSN<-distr_SkewNorm_EVSK_Th(omega,alpha)Compare the distinct elements of the estimated skewness vector and
the theoretical ones using indx_Elimination.
EsMSN$estSkew[indx_Elimination(3,3)]
#> [1] 0.8160645703 0.3872013315 0.0247111595 0.0919462860 -0.0532704409
#> [6] 0.0001422415 0.0991293326 0.0079415309 -0.0032085772 -0.0743442849
ThMSN$SkewX[indx_Elimination(3,3)]
#> [1] 0.68927167 0.34463583 0.00000000 0.17231792 0.00000000 0.00000000
#> [7] 0.08615896 0.00000000 0.00000000 0.00000000If one wishes to recover the estimated univariate skewness and
kurtosis of the components \(X1\),
\(X2\) and \(X3\) of \(X\), then, using
indx_UnivMomCum,
EsMSN$estSkew[indx_UnivMomCum(3,3)] ## Get univariate skewness for X1,X2,X3
#> [1] 0.81606457 0.09912933 -0.07434428
EsMSN$estKurt[indx_UnivMomCum(3,4)] ## Get univariate kurtosis for X1,X2,X3
#> [1] 1.2089545 -0.1196923 0.1323258An 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.
Esti_Skew_Mardia(data)
#> $Mardia.Skewness
#> [1] 1.175531
#>
#> $p.value
#> [1] 1.142991e-36
ThMSN$SkewX.tot
#> [1] 0.9279208The 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 […]
This work has been partially supported by the project TKP2021-NKTA of the University of Debrecen, Hungary.
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.