---
title: "MultiStatM: overview"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Overview}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  markdown: 
    wrap: 72
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(MultiStatM)
```

## 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")`.

+------------------------+-----------------------------+---------------+
| **MultiStatM 1.2.1**   | **MultiStatM 2.0.0**        | **Family**    |
+========================+=============================+===============+
|conv_Cum2Mom            | Cum2Mom,                    | Moments and   |
|conv_Cum2MomMulti       | Type=Univariate/Multivariate| cumulants     |
+------------------------+-----------------------------+---------------+
|conv_Mom2Cum            | Mom2Cum,                    | Moments and   |
|conv_Mom2CumMulti       | Type=Univariate/Multivariate| 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   | MomCumZabs,                 | Moments and   |
|distr_Zabs_MomCum_Th    | Type=Univariate/Multivariate| cumulants     |
+------------------------+-----------------------------+---------------+
|distr_SkewNorm_EVSK_Th  | EVSKSkewNorm                | Moments and   |
|                        |                             | cumulants     |
+------------------------+-----------------------------+---------------+
|distr_Uni_EVSK_Th       | EVSKUniS                    | Moments and   |
|distr_UniAbs_EVSK_Th    | Type=Standard/Modulus       | 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        | SampleKurt                  | Estimation    |
|Esti_Kurt_MRSz          | Type=Mardia/MRSz/Total      |               |
|Esti_Kurt_Total         |                             |               |
+------------------------+-----------------------------+---------------+
|Esti_Skew_Mardia        | SampleSkew                  | Estimation    |
| Esti_Skew_MRSz         | Type=Mardia/MRSz            |               |
+------------------------+-----------------------------+---------------+
|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           | HermiteCoeff                | Hermite       |
|Hermite_CoeffMulti      | Type=Univariate/Multivariate| Polynomials   |
+------------------------+-----------------------------+---------------+
|Hermite_Poly_HN         | HermiteN                    | Hermite       |
|Hermite_Poly_HN_Multi   | Type=Univariate/Multivariate| Polynomials   |
+------------------------+-----------------------------+---------------+
|Hermite_Poly_NH_Inv     | HermiteN2X                  | Hermite       |
|Hermite_Poly_NH_Multi_In| Type=Univariate/Multivariate| Polynomials   |
+------------------------+-----------------------------+---------------+
| Hermite_Nth            | Eliminated: use HermiteN    |               |
+------------------------+-----------------------------+---------------+
| Hermite_N_Cov_X1_X2    | HermiteCov12                | Hermite       |
|                        |                             | Polynomials   |
+------------------------+-----------------------------+---------------+
|indx_Commutator_Kmn     | CommutatorIndx              | Commutators   |
|indx_Commutator_Kperm   | Type=Kmn/Kperm/Mixing/Moment|               |
|indx_Commutator_Mixing  |                             |               |
|indx_Commutator_Moment  |                             |               |
+------------------------+-----------------------------+---------------+
| indx_Elimination       | EliminIndx                  | Commutators   |
+------------------------+-----------------------------+---------------+
| indx_Qplication        | QplicIndx                   | Commutators   |
+------------------------+-----------------------------+---------------+
| indx_Symmetry          | SymIndx                     | Commutators   |
+------------------------+-----------------------------+---------------+
| indx_UnivMomCum        | UnivMomCum                  | Commutators   |
+------------------------+-----------------------------+---------------+
| matr_Commutator_Kmn    | CommutatorMatr              | Commutators   |
| matr_Commutator_Kperm  | Type=Kmn/Kperm/Mixing/Moment|               |
| matr_Commutator_Mixing |                             |               |
| matr_Commutator_Moment |                             |               |
+------------------------+-----------------------------+---------------+
| matr_Elimination       | EliminMatr                  | Commutators   |
+------------------------+-----------------------------+---------------+
| matr_Qplication        | QplicMatr                   | Commutators   |
+------------------------+-----------------------------+---------------+
| matr_Symmetry          | SymMatr                     | Commutators   |
+------------------------+-----------------------------+---------------+
|Partition_2Perm         | Partitions                  | Partitions    |
|Partition_Diagrams      | Type=2Perm/Diagram/Indecomp |               |              |        ClosedNoLoops   |      /Pairs                 |               |              
|Partition_Indecomposable|                             |               |
|Partition_Pairs         |                             |               |
+------------------------+-----------------------------+---------------+
|Permutation_Inverse     | PermutationInv              | Partitions    |
+------------------------+-----------------------------+---------------+
|Partition_Type_All      | PartitionTypeAll            | Partitions    |
+------------------------+-----------------------------+---------------+



## 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}$\texttt{S\_N\_r}$\left[ r\right] +1$ up to $\sum_{r\leq R}$\texttt{S\_N\_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

```{r}
PTA<-PartitionTypeAll(4)
```

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

```{r}
PTA$S_N_r
```

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

```{r}
PTA$Part.class[[2]]
```

The second partition matrix of size 3 is

```{r}
PTA$Part.class[[10]] 
```

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

```{r}
PTA$S_r_j[[2]] 
```

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`:

```{r}
PTA$eL_r
```

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$.


```{r}
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) 
a3%x%a1%x%a2 
```

The same result can be obtained by using `CommutatorIndx`


```{r}
p1[CommutatorIndx(Type="Kperm",c(3,1,2),c(2,3,2))]
```

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.

```{r}
a<-c(1,2)
a3<-a%x%a%x%a
a3
c(EliminMatr(2,3)%*%a3)
c(QplicMatr(2,3)%*%EliminMatr(2,3)%*%a3)
```

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


```{r}
x<-c(1,2,3)
H2<-HermiteN(x,N=2,Type="Multivariate")
H2[[1]]
H2[[2]]
```


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


```{r}
H2<-HermiteN(x,Sig2=4*diag(3),N=2,Type="Multivariate")
H2[[1]]
H2[[2]]
```

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

```{r}
HermiteN2X(H2,N=2,Sig2=4*diag(3),Type="Multivariate")[[1]] 
```

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}}$.

```{r}
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

```{r}
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
```

Getting back to moments

```{r}
Cum2Mom(cum,Type="Multivariate")
```

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

```{r}
mu[[3]][EliminIndx(2,3)]
```

Alternatively one can also use the elimination matrix

```{r}
r.mu<-EliminMatr(2,3)%*% mu[[3]]
c(r.mu)
```

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.

```{r}
c(QplicMatr(2,3)%*%r.mu)
```

The same result by using `QplicIndx`

```{r}
r.mu[QplicIndx(2,3)]
```

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 \CRANpkg{MultiStatM}  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

```{r}
alpha<-c(10,5,0)
omega<-diag(3)
MSN<-MomCumSkewNorm(r=3,omega,alpha,nMu=TRUE)
round(MSN$Mu[[3]],3)
round(MSN$CumX[[3]],3)
```

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

```{r}
EVSKUniS(3,  Type="Standard")$Kurt.U
```

## 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.

```{r}
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`.

```{r}
EsMSN$estSkew[EliminIndx(3,3)]   
ThMSN$SkewX[EliminIndx(3,3)]   

```

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

```{r}
EsMSN$estSkew[UnivMomCum(3,3)]  ## Get univariate skewness for X1,X2,X3
EsMSN$estKurt[UnivMomCum(3,4)]  ## Get univariate kurtosis for X1,X2,X3
```

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`.

```{r}
SampleSkew(data,Type="Mardia")
ThMSN$SkewX.tot
```

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 [...]


```{r}
SampleSkew(data,Type="MRSz")
as.vector(t(c(diag(3))%x%diag(3))%*%ThMSN$SkewX)
```

```{r}
c(t(c(diag(3))%x%diag(3))%*%ThMSN$SkewX)  ## Theoretical MRS skewness vector
```

## 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.