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
skewnormal and the skewt 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 higherorder 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 wellknown 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 k1}\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 higherorder cumulants leads to highdimensional 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 (OuldBaba (2015), which goes as far as the sixthorder 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")
.
MultiStatM 1.2.1conv_Cum2Mom conv_Cum2MomMulti 
MultiStatM 2.0.0Cum2Mom, Type=Univariate/Multivariate 
FamilyMoments and cumulants 

conv_Mom2Cum conv_Mom2CumMulti  Mom2Cum, Type=Univariate/Multivariate  Moments and cumulants  
conv_Stand_Multi  MVStandardize  
distr_CFUSN_MomCum_Th  MomCumCFUSN  Moments and cumulants  
distr_SkewNorm_MomCum_Th  MomCumSkewNorm  Moments and cumulants  
distr_Uni_MomCum_Th  MomCumUniS  Moments and cumulants  
distr_ZabsM_MomCum_Th distr_Zabs_MomCum_Th  MomCumZabs, Type=Univariate/Multivariate  Moments and cumulants  
distr_SkewNorm_EVSK_Th  EVSKSkewNorm  Moments and cumulants  
distr_Uni_EVSK_Th distr_UniAbs_EVSK_Th  EVSKUniS Type=Standard/Modulus  Moments and cumulants  
distr_CFUSN_Rand  rCFUSN,  Random Generation  
distr_CFUSSD_Rand  rCFUSSD  Random Generation  
distr_SkewNorm_Rand  rSkewNorm  Random Generation  
distr_Uni_Rand  rUniS  Random Generation  
Esti_Kurt_Mardia Esti_Kurt_MRSz Esti_Kurt_Total  SampleKurt Type=Mardia/MRSz/Total  Estimation  
Esti_Skew_Mardia Esti_Skew_MRSz  SampleSkew Type=Mardia/MRSz  Estimation  
Esti_Kurt_Variance_Th  VarianceKurt  Estimation  
Esti_Skew_Variance_Th  VarianceSkew  Estimation  
Esti_EVSK  SampleEVSK  Estimation  
Esti_Hermite_Poly_HN_M  SampleHermiteN  Estimation   
Esti_Gram_Charlier  SampleGC  Estimation  
Esti_MMom_MCum  SampleMomCum  Estimation  
Esti_Variance_Skew_Kurt  SampleVarianceSkewKurt  
Hermite_Coeff Hermite_CoeffMulti  HermiteCoeff Type=Univariate/Multivariate  Hermite Polynomials  
Hermite_Poly_HN Hermite_Poly_HN_Multi  HermiteN Type=Univariate/Multivariate  Hermite Polynomials  
Hermite_Poly_NH_Inv Hermite_Poly_NH_Multi_In  HermiteN2X Type=Univariate/Multivariate  Hermite Polynomials  
Hermite_Nth  Eliminated: use HermiteN  
Hermite_N_Cov_X1_X2  HermiteCov12  Hermite Polynomials  
indx_Commutator_Kmn indx_Commutator_Kperm indx_Commutator_Mixing indx_Commutator_Moment  CommutatorIndx Type=Kmn/Kperm/Mixing/Moment  Commutators  
indx_Elimination  EliminIndx  Commutators  
indx_Qplication  QplicIndx  Commutators  
indx_Symmetry  SymIndx  Commutators  
indx_UnivMomCum  UnivMomCum  Commutators  
matr_Commutator_Kmn matr_Commutator_Kperm matr_Commutator_Mixing matr_Commutator_Moment  CommutatorMatr Type=Kmn/Kperm/Mixing/Moment  Commutators  
matr_Elimination  EliminMatr  Commutators  
matr_Qplication  QplicMatr  Commutators  
matr_Symmetry  SymMatr  Commutators  
Partition_2Perm Partition_Diagrams Partition_Indecomposable Partition_Pairs  Partitions Type=2Perm/Diagram/Indecomp  Partitions  
Permutation_Inverse  PermutationInv  Partitions  
Partition_Type_All  PartitionTypeAll  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 nonempty 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}}_{\{rl\}}\). 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 columnsums \(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
S_N_r
provides the number of partitions with
r=1
to r=4
blocks:
All the partition matrices are listed in Part.class
, for
example the first partition matrix of size 2 among 7 is
The second partition matrix of size 3 is
etc.. If one interested in the number of partitions with different
types for r=2
then consider the list S_r_j
,
i.e.
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.
The CommutatorMatr
function produces commutators and
selection matrices. The use of matrices allows represent as linear
combinations problems of permutation and powers of Tproducts. 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
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 qplicated elements in powers of Tproducts.
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 dvector X from available vector of Tmoments and
Tcumulants. The latter eliminates the duplicated/qplicated elements in
a Tvector 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 Tproduct 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 Tproduct 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.
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 Tderivative 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 THermite 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}_{n1}(\mathbf{X}) \otimes \mathbf{X} (n1) \mathbf{H}_{n2}(\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:
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}}\).
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 Tderivative 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 viceversa 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:nr+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)^{r1} (r1)!\sum_{\sum l_j =r, \sum j l_j = n} \prod_{j=1:nr+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 2variate
standard normal distribution with null mean vector and covariance matrix
with unit elements on the main diagonal and offdiagonal 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
Alternatively one can also use the elimination matrix
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.
The same result by using QplicIndx
The MomCum
functions provide theoretical moments and
cumulants for some common multivariate distributions: Skewnormal,
Canonical Fundamental Skewnormal (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}_{d1}\).
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
variancecovariance 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 skewnormal distribution introduced by @azzalini1996multivariate, whose marginal
densities are scalar skewnormals. A \(d\)dimensional random vector \(\mathbf{X}\) is said to have a multivariate
skewnormal 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
rth 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 skewnormal 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. ArellanoValle 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 skewnormal 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
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 thirdorder Tcumulant (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 thirdorder central moments and thirdorder 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 perrow, 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 fourthorder Tcumulant 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 dataestimated 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 skewnormal
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.
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
pvalue 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.345182e23
ThMSN$SkewX.tot
#> [1] 0.9279208
The MRS skewness vector and index are provided together with the pvalue 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 TKP2021NKTA of the University of Debrecen, Hungary.
ArellanoValle, 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 SkewNormal and Related Distributions such as the Skewt and the SUN (version 2.0.2). URL http://azzalini.stat.unipd.it/SN/, https://cran.rproject.org/package=sn
Azzalini, A,, Dalla Valle, A. (1996) The multivariate skewnormal 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), 959974.
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.Rproject.org/package=kStatistics
Franceschini, C. and Loperfido, N. (2017a). MultiSkew: Measures, Tests and Removes Multivariate Skewness. R package version 1.1.1. https://CRAN.Rproject.org/package=MultiSkew
Franceschini, C. and Loperfido, N. (2017b). MaxSkew: Orthogonal Data Projections with Maximal Skewness. R package version 1.1. https://CRAN.Rproject.org/package=MaxSkew
Holmquist, B. (1996). The dvariate 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), 607644.
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), 708728.
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.Rproject.org/package=moments
Novomestky, F. (2021). matrixcalc: Collection of Functions for Matrix Calculations. R package version 1.05. https://CRAN.Rproject.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.
OuldBaba, 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.