Multivariate Geometric Skew-Normal Distribution

TL;DR

This paper introduces the multivariate geometric skew-normal distribution, explores its properties, and demonstrates its advantages and ease of use through simulations and real data analysis.

Contribution

It extends the geometric skew-normal distribution to multivariate cases, providing new properties, characterization results, and an EM algorithm for parameter estimation.

Findings

The distribution can take various shapes in its joint density.

The EM algorithm performs well in parameter estimation.

The model fits real data effectively.

Abstract

Azzalini (1985) introduced a skew-normal distribution of which normal distribution is a special case. Recently Kundu (2014) introduced a geometric skew-normal distribution and showed that it has certain advantages over Azzalini's skew-normal distribution. In this paper we discuss about the multivariate geometric skew-normal distribution. It can be used as an alternative to Azzalini's skew normal distribution. We discuss different properties of the proposed distribution. It is observed that the joint probability density function of the multivariate geometric skew normal distribution can take variety of shapes. Several characterization results have been established. Generation from a multivariate geometric skew normal distribution is quite simple, hence the simulation experiments can be performed quite easily. The maximum likelihood estimators of the unknown parameters can be obtained quite conveniently using expectation maximization (EM) algorithm. We perform some simulation experiments and it is observed that the performances of the proposed EM algorithm are quite satisfactory. Further, the analyses of two data sets have been performed, and it is observed that the proposed methods and the model work very well.

Tables

Table 1. Table 1: Average estimates and MSEs of 𝝁 ^ ^ 𝝁 \widehat{\mbox{\boldmath$\mu$}} and 𝚺 ^ ^ 𝚺 \widehat{\mbox{\boldmath$\Sigma$}} , when p 𝑝 p = 0.5 and it is known

	0.0000	0.0000	1.0000	1.0000
$𝝁$	0.0053	0.0047	1.0097	1.0055
	(0.0984)	(0.1213)	(0.1430)	(0.1237)
	2.0000	2.0000	1.0000	0.0000
	2.0024	1.9984	0.9942	-0.0060
	(0.3299)	(0.3602)	(0.2879)	(0.2262)
	2.0000	3.0000	2.0000	1.0000
	1.9984	2.9922	1.9907	0.9907
${\sigma }_{ij}$	(0.3602)	(0.4932)	(0.4130)	(0.3092)
	1.0000	2.0000	3.0000	2.0000
	0.9942	1.9907	2.9757	1.9823
	(0.2879)	(0.4130)	(0.5350)	(0.4054)
	0.0000	1.0000	2.0000	2.0000
	-0.0060	0.9907	1.9823	1.9859
	(0.2262)	(0.3092)	(0.4054)	(0.3675)

Table 2. Table 2: Average estimates and MSEs of 𝝁 ^ ^ 𝝁 \widehat{\mbox{\boldmath$\mu$}} , 𝚺 ^ ^ 𝚺 \widehat{\mbox{\boldmath$\Sigma$}} and p ^ ^ 𝑝 \widehat{p} , when p 𝑝 p = 0.5 and it is unknown

	0.0000	0.0000	1.0000	1.0000
$𝝁$	-0.0067	-0.0079	1.0094	1.0122
	(0.1029)	(0.1255)	(0.1553)	(0.1390)
	2.0000	2.0000	1.0000	0.0000
	2.0105	2.0103	1.0055	0.0011
	(0.3579)	(0.3964)	(0.3165)	(0.2377)
	2.0000	3.0000	2.0000	1.0000
	2.0103	3.0194	2.0095	1.0091
${\sigma }_{ij}$	(0.3964)	(0.5420)	(0.4443)	(0.3233)
	1.0000	2.0000	3.0000	2.0000
	1.0055	2.0095	2.9987	2.0006
	(0.3165)	(0.4443)	(0.5577)	(0.4161)
	0.0000	1.0000	2.0000	2.0000
	0.0011	1.0091	2.0006	2.0058
	(0.2377)	(0.3233)	(0.4161)	(0.3827)
$p$	0.5000
	0.5068
	(0.0433)

Table 3. Table 3: Average estimates and MSEs of 𝝁 ^ ^ 𝝁 \widehat{\mbox{\boldmath$\mu$}} and 𝚺 ^ ^ 𝚺 \widehat{\mbox{\boldmath$\Sigma$}} , when p 𝑝 p = 0.75 and it is known

	0.0000	0.0000	1.0000	1.0000
$𝝁$	0.0057	0.0046	1.0118	1.0068
	(0.1205)	(0.1481)	(0.1605)	(0.1350)
	2.0000	2.0000	1.0000	0.0000
	2.0015	1.9969	0.9923	-0.0065
	(0.3126)	(0.3416)	(0.2720)	(0.2099)
	2.0000	3.0000	2.0000	1.0000
	1.9969	2.9906	1.9887	0.9898
${\sigma }_{ij}$	(0.3416)	(0.4645)	(0.3875)	(0.2907)
	1.0000	2.0000	3.0000	2.0000
	0.9923	1.9897	2.9778	1.9864
	(0.2720)	(0.3875)	(0.4848)	(0.3864)
	0.0000	1.0000	2.0000	2.0000
	-0.0065	0.9898	1.9864	1.9903
	(0.2099)	(0.2907)	(0.3684)	(0.3320)

Table 4. Table 4: Average estimates and MSEs of 𝝁 ^ ^ 𝝁 \widehat{\mbox{\boldmath$\mu$}} , 𝚺 ^ ^ 𝚺 \widehat{\mbox{\boldmath$\Sigma$}} and p ^ ^ 𝑝 \widehat{p} , when p 𝑝 p = 0.75 and it is unknown

	0.0000	0.0000	1.0000	1.0000
$𝝁$	0.0047	0.0035	1.0141	1.0093
	(0.1392)	(0.1704)	(0.1938)	(0.1615)
	2.0000	2.0000	1.0000	0.0000
	2.0183	2.0258	1.0263	0.0121
	(0.3793)	(0.4260)	(0.3453)	(0.2620)
	2.0000	3.0000	2.0000	1.0000
	2.0258	3.0351	2.0272	1.0088
	(0.4260)	(0.5807)	(0.4842)	(0.3448)
	1.0000	2.0000	3.0000	2.0000
${\sigma }_{ij}$	1.0263	2.0272	3.0120	1.9961
	(0.3453)	(0.4842)	(0.5897)	(0.4256)
	0.0000	1.0000	2.0000	2.0000
	0.0121	1.0088	1.9961	1.9892
	(0.2620)	(0.3448)	(0.4256)	(0.3726)
$p$	0.7500
	0.7575
	(0.0446)

Table 5. Table 5: Four different stiffness measurements of 30 boards

No.	$x_1$	$x_2$	$x_3$	$x_4$	No.	$x_1$	$x_2$	$x_3$	$x_4$
1	1889	1651	1561	1778	2	2403	2048	2087	2197
3	2119	1700	1815	2222	4	1645	1627	1110	1533
5	1976	1916	1614	1883	6	1712	1712	1439	1546
7	1943	1685	1271	1671	8	2104	1820	1717	1874
9	2983	2794	2412	2581	10	1745	1600	1348	1508
11	1710	1591	1518	1667	12	2046	1907	1627	1898
13	1840	1841	1595	1741	14	1867	1685	1493	1678
15	1859	1649	1389	1714	16	1954	2149	1180	1281
17	1325	1170	1002	1176	18	1419	1371	1251	1308
19	1828	1634	1602	1755	20	1725	1594	1313	1646
21	2276	2189	1547	2111	22	1899	1614	1422	1477
23	1633	1513	1290	1516	24	2061	1867	1646	2037
25	1856	1493	1356	1533	26	1727	1412	1238	1469
27	2168	1896	1701	1834	28	1655	1675	1414	1597
29	2326	2301	2065	2234	30	1490	1382	1214	1284

Equations

$$f(x;\mu ,\sigma ,\lambda )=\frac{2}{\sigma }\varphi \left(\frac{x-\mu }{\sigma }\right)\Phi \left(\frac{\lambda (x-\mu )}{\sigma }\right),-\infty <x,\mu ,\lambda <\infty ,\sigma >0,$$

$$g({\mbox{\boldmath$z$}})=2\phi_{d}({\mbox{\boldmath$z$}};{\mbox{\boldmath$\Omega$}})\Phi({\mbox{\boldmath$\alpha$}}^{T}{\mbox{\boldmath$z$}}),\ \ \ {\mbox{\boldmath$z$}}\in\mathbb{R}^{d},$$

$$X\stackrel{dist}{=}\sum _{i=1}^N{X}_i$$

$$f_X(x;\mu ,\sigma ,p)=\sum _{k=1}^{\infty }\frac{p}{\sigma \sqrt{k}}\varphi \left(\frac{x-k\mu }{\sigma \sqrt{k}}\right)(1-p{)}^{k-1}.$$

$$M_X(t)=\frac{p{e}^{\mu t+\frac{{\sigma }^2{t}^2}{2}}}{1-(1-p)e^{\mu t+\frac{{\sigma }^2{t}^2}{2}}},t\in A_1(\mu ,\sigma ,p),$$

$$A_1(\mu ,\sigma ,p)$$

$$K_X(t)=\ln M_X(t)=\ln p+\mu t+\frac{{\sigma }^2{t}^2}{2}-\ln \left(1-(1-p)e^{\mu t+\frac{{\sigma }^2{t}^2}{2}}\right).$$

$$E(X)$$

$$V(X)$$

$${\gamma }_1$$

$${\gamma }_2$$

$${\mbox{\boldmath$X$}}\stackrel{{\scriptstyle dist}}{{=}}\sum_{i=1}^{N}{\mbox{\boldmath$X$}}_{i},$$

$$F_{\mbox{\boldmath$X$}}({\mbox{\boldmath$x$}};{\mbox{\boldmath$\mu$}},{\mbox{\boldmath$\Sigma$}},p)=\sum_{k=1}^{\infty}p(1-p)^{k-1}\Phi_{d}({\mbox{\boldmath$x$}};k{\mbox{\boldmath$\mu$}},k{\mbox{\boldmath$\Sigma$}})$$

$$\displaystyle f_{\mbox{\boldmath$X$}}({\mbox{\boldmath$x$}};{\mbox{\boldmath$\mu$}},{\mbox{\boldmath$\Sigma$}},p)$$

$$M_{\mbox{\boldmath$X$}}({\mbox{\boldmath$t$}})=\frac{pe^{{\mbox{\boldmath$\mu$}}^{T}{\mbox{\boldmath$t$}}+\frac{1}{2}{\mbox{\boldmath$t$}}^{T}{\mbox{\boldmath$\Sigma$}}{\mbox{\boldmath$t$}}}}{1-(1-p)e^{{\mbox{\boldmath$\mu$}}^{T}{\mbox{\boldmath$t$}}+\frac{1}{2}{\mbox{\boldmath$t$}}^{T}{\mbox{\boldmath$\Sigma$}}{\mbox{\boldmath$t$}}}},\ \ \ \ {\mbox{\boldmath$t$}}\in A_{d}({\mbox{\boldmath$\mu$}},{\mbox{\boldmath$\Sigma$}},p),$$

$$\displaystyle A_{d}({\mbox{\boldmath$\mu$}},{\mbox{\boldmath$\Sigma$}},p)$$

$${\mbox{\boldmath$X$}}=\left(\matrix{{\mbox{\boldmath$X$}}_{1}\cr{\mbox{\boldmath$X$}}_{2}\cr}\right),\ \ {\mbox{\boldmath$\mu$}}=\left(\matrix{{\mbox{\boldmath$\mu$}}_{1}\cr{\mbox{\boldmath$\mu$}}_{2}\cr}\right),\ \ \ {\mbox{\boldmath$\Sigma$}}=\left(\matrix{{\mbox{\boldmath$\Sigma$}}_{11}&{\mbox{\boldmath$\Sigma$}}_{12}\cr{\mbox{\boldmath$\Sigma$}}_{21}&{\mbox{\boldmath$\Sigma$}}_{22}\cr}\right).$$

$${\mbox{\boldmath$X$}}_{2}\sim MGSN_{d-h}(p,{\mbox{\boldmath$\mu$}}_{2},{\mbox{\boldmath$\Sigma$}}_{22}).$$

$$\displaystyle M_{\mbox{\boldmath$Z$}}({\mbox{\boldmath$t$}})$$

$$A_s^D$$

$$\left.E(X_{i})=\frac{\partial}{\partial t_{i}}M_{\mbox{\boldmath$X$}}({\mbox{\boldmath$t$}})\right|_{{\mbox{\boldmath$t$}}={\mbox{\boldmath$0$}}}=\frac{\mu_{i}}{p}$$

$$\left.E(X_{i}X_{j})=\frac{\partial^{2}}{\partial t_{i}\partial t_{j}}M_{\mbox{\boldmath$X$}}({\mbox{\boldmath$t$}})\right|_{{\mbox{\boldmath$t$}}={\mbox{\boldmath$0$}}}=\frac{p\sigma_{ij}+\mu_{i}\mu_{j}(2-p)}{p^{2}}.$$

$$Var(X_i)=\frac{p{\sigma }_{ii}+{\mu }_i^2(1-p)}{p^2},$$

$$Cov(X_i,X_j)=\frac{p{\sigma }_{ij}+{\mu }_i{\mu }_j(1-p)}{p^2},$$

$$Corr(X_i,X_j)=\frac{p{\sigma }_{ij}+{\mu }_i{\mu }_j(1-p)}{\sqrt{p{\sigma }_{ii}+{\mu }_i^2(1-p)}\sqrt{p{\sigma }_{jj}+{\mu }_j^2(1-p)}}.$$

$${\mu }_{i_1,\dots ,i_s}^{(r_1,\dots ,r_s)}=E\left[\prod _{k=1}^s(X_{r_k}-{\mu }_{r_k}{)}^{i_k}\right],$$

$${\beta }_1=\sum _{r,s,t=1}^d\sum _{r^{\prime },s^{\prime },t^{\prime }=1}^d{\sigma }^{r{r}^{\prime }}{\sigma }^{s{s}^{\prime }}{\sigma }^{t{t}^{\prime }}{\mu }_{111}^{rst}{\mu }_{111}^{r^{\prime }{s}^{\prime }{t}^{\prime }},$$

$${\mu }_{111}^{lhm}=\frac{1}{p^4}\left\{p(1-p)(2-p){\mu }_h{\mu }_l{\mu }_m+p^2(1-p)({\mu }_m{\sigma }^{hl}+{\mu }_l{\sigma }^{hm}+{\mu }_h{\sigma }^{lm})\right\}.$$

$$p\ {\mbox{\boldmath$\mu$}}_{\mbox{\boldmath$X$}}={\mbox{\boldmath$\mu$}}\ \ \ \ \hbox{and}\ \ \ \ p^{2}\ {\mbox{\boldmath$\Sigma$}}_{\mbox{\boldmath$X$}}=p\ {\mbox{\boldmath$\Sigma$}}+(1-p)\ {\mbox{\boldmath$\mu$}}{\mbox{\boldmath$\mu$}}^{T}.$$

$$\left|\matrix{-\lambda{\mbox{\boldmath$\Sigma$}}_{11}&{\mbox{\boldmath$\Sigma$}}_{12}\cr{\mbox{\boldmath$\Sigma$}}_{21}&-\lambda{\Sigma}_{22}\cr}\right|=0.$$

$$\left(\matrix{{\mbox{\boldmath$\alpha$}}^{T}{\mbox{\boldmath$X$}}_{1}\cr{\mbox{\boldmath$\beta$}}^{T}{\mbox{\boldmath$X$}}_{2}\cr}\right)\sim\hbox{MGSN}\left(2,p,\left(\matrix{{\mbox{\boldmath$\alpha$}}^{T}{\mbox{\boldmath$\mu$}}_{1}\cr{\mbox{\boldmath$\beta$}}^{T}{\mbox{\boldmath$\mu$}}_{2}\cr}\right),\left(\matrix{{\mbox{\boldmath$\alpha$}}^{T}{\mbox{\boldmath$\Sigma$}}_{11}{\mbox{\boldmath$\alpha$}}&{\mbox{\boldmath$\alpha$}}^{T}{\mbox{\boldmath$\Sigma$}}_{12}{\mbox{\boldmath$\beta$}}\cr{\mbox{\boldmath$\beta$}}^{T}{\mbox{\boldmath$\Sigma$}}_{21}{\mbox{\boldmath$\alpha$}}&{\mbox{\boldmath$\beta$}}^{T}{\mbox{\boldmath$\Sigma$}}_{22}{\mbox{\boldmath$\beta$}}\cr}\right)\right).$$

$$\hbox{corr}({\mbox{\boldmath$\alpha$}}^{T}{\mbox{\boldmath$X$}}_{1},{\mbox{\boldmath$\beta$}}^{T}{\mbox{\boldmath$X$}}_{2})=\frac{p{\mbox{\boldmath$\alpha$}}^{T}{\mbox{\boldmath$\Sigma$}}_{12}{\mbox{\boldmath$\beta$}}^{T}}{\sqrt{p{\mbox{\boldmath$\alpha$}}^{T}{\mbox{\boldmath$\Sigma$}}_{11}{\mbox{\boldmath$\alpha$}}}\sqrt{p{\mbox{\boldmath$\beta$}}^{T}{\mbox{\boldmath$\Sigma$}}_{22}{\mbox{\boldmath$\beta$}}}}={\mbox{\boldmath$\alpha$}}^{T}{\mbox{\boldmath$\Sigma$}}_{12}{\mbox{\boldmath$\beta$}}^{T},$$

Full text

Key Words and Phrases: Skew-normal distribution; moment generating function; infinite divisible distribution; maximum likelihood estimators; EM algorithm; Fisher information matrix.

AMS 2000 Subject Classification: Primary 62F10; Secondary: 62H10

1Department of Mathematics and Statistics, Indian Institute of Technology Kanpur, Kanpur, Pin 208016, India.   e-mail: kundu@iitk.ac.in.

1 Introduction

Azzalini [3] proposed a class of three-parameter skew-normal distributions which includes the normal one. Azzalini’s skew normal (ASN) distribution has received a considerable attention in the last two decades due to its flexibility and its applications in different fields. The probability density function (PDF) of ASN takes the following form:

[TABLE]

where $\phi(x)$ and $\Phi(x)$ denote the standard normal PDF and standard normal cumulative distribution function (CDF), respectively, at the point $x$. Here $\mu$, $\sigma$ and $\lambda$ are known as the location, scale and skewness or tilt parameters, respectively. ASN distribution has an unimodal PDF, and it can be both positively or negatively skewed depending on the skewness parameter. Arnold and Beaver [2] provided an interesting interpretation of this model in terms of hidden truncation. This model has been used quite effectively to analyze skewed data in different fields due to its flexibility.

Later Azzalini and Dalla Valle [5] constructed a multivariate distribution with skew normal marginals. From now on we call it as Azzalini’s multivariate skew-normal (AMSN) distribution, and it can be defined as follows. A random vector {\mbox{\boldmathZ}}=(Z_{1},\ldots,Z_{d})^{T} is a $d$-dimensional AMSN distribution, if it has the following PDF

[TABLE]

where \phi_{d}({\mbox{\boldmathz}},{\mbox{\boldmath\Omega}}) denotes the PDF of the $d$-dimensional multivariate normal distribution with standardized marginals, and correlation matrix $\Omega$. We denote such a random vector as {\mbox{\boldmathZ}}\sim SN{}_{d}({\mbox{\boldmath\Omega}},{\mbox{\boldmath\alpha}}). Here the vector $\alpha$ is known as the shape vector, and it can be easily seen that the PDF of AMSN distribution is unimodal and can take different shapes. It has several interesting properties, and it has been used quite successfully to analyze several multivariate data sets in different areas because of its flexibility.

Although ASN distribution is a very flexible distribution, it cannot be used to model moderate or heavy tail data; see for example Azzalini and Capitanio [4]. It is well known to be a thin tail distribution. Since the marginals of AMSN are ASN, multivariate heavy tail data cannot be modeled by using AMSN. Due to this reason several other skewed distributions, often called skew-symmetric distributions, have been suggested in the literature using different kernel functions other than the normal kernel function and using the same technique as Azzalini [3]. Depending on the kernel function the resulting distribution can have moderate or heavy tail behavior. Among different such distributions, skew-t distribution is quite commonly used in practice, which can produce heavy tail distribution depending on the degrees of freedom of the associated $t$-distribution. It has a multivariate extension also. For a detailed discussions on different skew-symmetric distribution, the readers are referred to the excellent monograph by Azzalini and Capitanio [4].

Although ASN model is a very flexible one dimensional model, and it has several interesting properties, it is well known that computing the maximum likelihood estimators (MLEs) of the unknown parameters of an ASN model is a challenging issue. Azzalini [3] has shown that there is a positive probability that the MLEs of the unknown parameters of a ASN model do not exist. If all the data points have same sign, then the MLEs of unknown parameters of the ASN model may not exist. The problem becomes more severe for AMSN model, and the problem exists for other kernels also.

Recently, the author [9] proposed a new three-parameter skewed distribution, of which normal distribution is a special case. The proposed distribution can be obtained as a geometric sum of independent identically distributed (i.i.d.) normal random variables, and it is called as the geometric skew normal (GSN) distribution. It can be used quite effectively as an alternative to an ASN distribution. It is observed that the GSN distribution is a very flexible distribution, as its PDF can take different shapes depending on the parameter values. Moreover, the MLEs of the unknown parameters can be computed quite conveniently using the EM algorithm. It can be easily shown that the ‘pseudo-log-likelihood’ function has a unique maximum, and it can be obtained in explicit forms. Several interesting properties of the GSN distribution have also been developed by Kundu [9].

The main aim of this paper is to consider the multivariate geometric skew-normal (MGSN) distribution, develop its various properties and discuss different inferential issues. Several characterization results and dependence properties have also been established. It is observed that the generation from a MGSN distribution is quite simple, hence simulation experiments can be performed quite conveniently. Note that the $d$-dimensional MGSN model has $d+1+d(d+1)/2$ unknown parameters. The MLEs of the unknown parameters can be obtained by solving $d+1+d(d+1)/2$ non-linear equations. We propose to use EM algorithm, and it is observed that the ’pseudo-log-likelihood’ function has a unique maximum, and it can be obtained in explicit forms. Hence, the implementation of the EM algorithm is quite simple, and the algorithm is very efficient. We perform some simulation experiments to see the performances of the proposed EM algorithm and the performances are quite satisfactory. We also perform the analyses of two data sets to illustrate how the proposed methods can be used in practice. It is observed that the proposed methods and the model work quite satisfactorily.

The main motivation to introduce the MGSN distribution can be stated as follows. Although there are several skewed distributions available in one-dimension, the same is not true in $\mathbb{R}^{d}$. The proposed MGSN distribution is a very flexible multivariate distribution which can produce variety of shapes. The joint PDF can be unimodal or multimodal and the marginals can have heavy tails depending on the parameters. It has several interesting statistical properties. Computation of the MLEs can be performed in a very simple manner even in high dimension. Hence, if it is known that the data are obtained from a multivariate skewed distribution, the proposed model can be used for analysis purposes. Generating random samples from a MGSN distribution is quite simple, hence any simulation experiment related to this distribution can be performed quite conveniently. Further, it is observed that in one of our data example the MLEs of AMSN do not exist, whereas the MLEs of MGSN distribution exist. Hence, in certain cases the implementation of MGSN distribution becomes easier than the AMSN distribution. The proposed MGSN distribution provides a choice to a practitioner of a new multivariate skewed distribution to analyze multivariate skewed data.

Rest of the paper is organized as follows. In Section 2, first we briefly describe the univariate GSN model, and discuss some of its properties, and then we describe MGSN model. Different properties are discussed in Section 3. In Section 4, we discuss the implementation of the EM algorithm, and some testing of hypotheses problems. Simulation results are presented in Section 5. The analysis of two data sets are presented in Section 6, and finally we conclude the paper in Section 7.

2 GSN and MGSN Distributions

We use the following notations in this paper. A normal random variable with mean $\mu$ and variance $\sigma^{2}$ will be denoted by N$(\mu,\sigma^{2})$. A $d$-variate normal random variable with mean vector $\mu$ and dispersion matrix $\Sigma$ will be denoted by N{}_{d}({\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}). The corresponding PDF and CDF at the point $x$ will be denoted by \phi_{d}({\mbox{\boldmathx}};{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}) and \Phi_{d}({\mbox{\boldmathx}};{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}), respectively. A geometric random variable with parameter $p$ will be denoted by GE$(p)$, and it has the probability mass function (PMF): $p(1-p)^{n-1}$ for $n=1,2,\ldots.$

2.1 GSN Distribution

Suppose $N\sim$ GE$(p)$ and $\{X_{i};i=1,2,\ldots,\}$ are i.i.d. Gaussian random variables. It is assumed that $N$ and $X_{i}$’s are independently distributed. Then the random variable

[TABLE]

is known as GSN random variable and its distribution will be denoted by GSN$(\mu,\sigma,p)$. Here, ‘$\displaystyle\stackrel{{\scriptstyle dist}}{{=}}$’ means equal in distribution. The GSN distribution can be seen as one of the compound geometric distributions. The PDF of $X$ takes the following form:

[TABLE]

When $\mu$ = 0 and $\sigma$ = 1, we say that $X$ has a standard GSN distribution, and it will be denoted by GSN$(p)$.

The standard GSN is symmetric about 0, and unimodal, but the PDF of GSN$(\mu,\sigma,p)$ can take different shapes. It can be unimodal or multimodal depending on $\mu$, $\sigma$ and $p$ values. The hazard function is always an increasing function. If $X\sim$ GSN$(\mu,\sigma,p)$, then the moment generating function (MGF) of $X$ becomes

[TABLE]

where

[TABLE]

The corresponding cumulant generating (CGF) function of $X$ is

[TABLE]

From (2), the mean, variance, skewness and kurtosis can be easily obtained as

[TABLE]

respectively. It is clear from the expressions of (3) and (4) that as $p\rightarrow 0$, $|E(X)|$ and $V(X)$ diverge to $\infty$. It indicates that GSN model can be used to model heavy tail data. It has been shown that the GSN law is infinitely divisible, and an efficient EM algorithm has been suggested to compute the MLEs of the unknown parameters.

2.2 MGSN Distribution

A $d$-variate MGSN distribution can be defined as follows. Suppose $N\sim$ GE$(p)$, \{{\mbox{\boldmathX}}_{i};i=1,2\ldots\} are i.i.d. N{}_{d}({\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}) random vectors and all the random variables are independently distributed. Define

[TABLE]

then $X$ is said to have a $d$-variate geometric skew-normal distribution with parameters $p$, $\mu$ and $\Sigma$, and its distribution will be denoted by MGSN{}_{d}(p,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}). If {\mbox{\boldmathX}}\sim MGSN{}_{d}(p,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}), then the CDF and PDF of $X$ become

[TABLE]

and

[TABLE]

respectively. Here \displaystyle\Phi_{d}({\mbox{\boldmathx}};k{\mbox{\boldmath\mu}},k{\mbox{\boldmath\Sigma}}) and \displaystyle\phi_{d}({\mbox{\boldmathx}};k{\mbox{\boldmath\mu}},k{\mbox{\boldmath\Sigma}}) denote the CDF and PDF of a $d$-variate normal distribution, respectively, with the mean vector \displaystyle k{\mbox{\boldmath\mu}} and dispersion matrix \displaystyle k{\mbox{\boldmath\Sigma}}.

If {\mbox{\boldmath\mu}}={\mbox{\boldmath0}} and {\mbox{\boldmath\Sigma}}={\mbox{\boldmathI}}, we say that $X$ is a standard $d$-variate MGSN random variable, and its distribution will be denoted by MGSN${}_{d}(p)$. The PDF of MGSN${}_{d}(p)$ is symmetric and unimodal, for all values of $d$ and $p$, whereas the PDF of MGSN{}_{d}(p,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}) may not be symmetric, and it can be unimodal or multimodal depending on parameter values. The MGF of MGSN can be obtained in explicit form. If {\mbox{\boldmathX}}\sim MGSN{}_{d}(p,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}), then the MGF of $X$ is

[TABLE]

where

[TABLE]

Further the generation of MGSN distribution is very simple. The following algorithm can be used to generate samples from a MGSN random variable.

Algorithm 1:

• •

Step 1: Generate $n$ from a GE$(p)$

• •

Step 2: Generate {\mbox{\boldmathX}}\sim N{}_{d}(n{\mbox{\boldmath\mu}},n{\mbox{\boldmath\Sigma}}).

In Figure 1 we provide the joint PDF of a bivariate geometric skew normal distribution for different parameter values: (a) $p$ = 0.75, $\mu_{1}=\mu_{2}$ = 0, $\sigma_{1}^{2}=\sigma_{2}^{2}$ = 2, $\sigma_{12}=\sigma_{21}=0$, (b) $p$ = 0.5, $\mu_{1}=\mu_{2}$ = 2.0, $\sigma_{1}^{2}=\sigma_{2}^{2}$ = 1, $\sigma_{12}=\sigma_{21}=-0.5$, (c) $p$ = 0.15, $\mu_{1}$ = 2.0, $\mu_{2}$ = 1.0, $\sigma_{1}^{2}=\sigma_{2}^{2}$ = 1, $\sigma_{12}=\sigma_{21}=-0.5$, (d) $p$ = 0.15, $\mu_{1}$ = 0.5, $\mu_{2}$ = -2.5, $\sigma_{1}^{2}=\sigma_{2}^{2}$ = 1.0, $\sigma_{12}=\sigma_{21}=0.5$.

3 Properties

In this section we discuss different properties of a MGSN distribution. We use the following notations:

[TABLE]

Here the vectors $X$ and $\mu$ are of the order $d$ each, and the matrix {\mbox{\boldmath\Sigma}}_{11} is of the order $h\times h$. Rest of the quantities are defined, so that they are compatible. The following result provides the marginals of a MGSN distribution.

Result 1: If \displaystyle{\mbox{\boldmathX}}\sim\hbox{MGSN}_{d}(p,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}) and \displaystyle{\mbox{\boldmathX}}_{1}\sim\hbox{MGSN}_{h}(p,{\mbox{\boldmath\mu}}_{1},{\mbox{\boldmath\Sigma}}_{11}) then

[TABLE]

Proof: The result easily follows from the MGF of MGSN as provided in (6).

We further have the following results similar to the multivariate normal distribution. The result may be used for testing simultaneously a set of linear hypothesis on the parameter vector $\mu$ or it may have some independent interest also; see for example Rao [12].

Theorem 1: If {\mbox{\boldmathX}}\sim MGSN{}_{d}(p,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}), then {\mbox{\boldmathZ}}={\mbox{\boldmathD}}{\mbox{\boldmathX}}\sim MGSN{}_{s}(p,{\mbox{\boldmathD}}{\mbox{\boldmath\mu}},{\mbox{\boldmathD}}{\mbox{\boldmath\Sigma}}{\mbox{\boldmathD}}^{T}), where $D$ is a $s\times d$ matrix of rank $s\leq d$.

Proof: The MGF of the random vector $Z$ is

[TABLE]

where

[TABLE]

Hence the result follows.

If {\mbox{\boldmathX}}=(X_{1},\ldots,X_{d})^{T}\sim MGSN{}_{d}(p,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}), and if we denote {\mbox{\boldmath\mu}}^{T}=(\mu_{1},\ldots,\mu_{d}), {\mbox{\boldmath\Sigma}}=((\sigma_{ij})), then the moments and cumulants of $X$, for $i,j=1,2,\ldots,d$, can be obtained from the MGF as follows:

[TABLE]

[TABLE]

Hence,

[TABLE]

[TABLE]

and

[TABLE]

It is clear from (11) that the correlation between $X_{i}$ and $X_{j}$ for $i\neq j$, not only depends on $\sigma_{ij}$, but it also depends on $\mu_{i}$ and $\mu_{j}$. For fixed $p$, $\sigma_{ij}$, if $\mu_{j}\rightarrow\infty$ and $\mu_{i}\rightarrow\infty$, then $Corr(X_{i},X_{j})\rightarrow$ 1, and if $\mu_{j}\rightarrow\infty$ and $\mu_{i}\rightarrow-\infty$, then $Corr(X_{i},X_{j})\rightarrow$ -1. From (11) it also follows that if $X$ is a standard $d$-variate MGSN random variable, i.e. for $i\neq j$, $\mu_{i}=\mu_{j}=\sigma_{ij}$ = 0, hence $Corr(X_{i},X_{j})$ = 0. Therefore, in this case although $X_{i}$ and $X_{j}$ are uncorrelated, they are not independent.

Now we would like to compute the multivariate skewness indices of the MGSN distribution. Different multivariate skewness measures have been introduced in the literature. Among them the skewness index of Mardia [10, 11] is the most popular one. To define Mardia’s multivariate skewness index let us introduce the following notations of a random vector {\mbox{\boldmathX}}=(X_{1},\ldots,X_{d}).

[TABLE]

where $\displaystyle\mu_{r_{k}}=E(X_{r_{k}}),k=1,\ldots,s$. Mardia [10] defined the multivariate skewness index as

[TABLE]

here $\sigma^{jk}$ for $j,k=1,\ldots,d$ denotes the $(j,k)$-th element of the inverse of the dispersion matrix of the random vector $X$. In case of MGSN distribution

[TABLE]

It is clear from (12) that if $p$ = 1 then $\beta_{1}$ = 0. Also if $\mu_{j}$ = 0 for all $j=1,\ldots,d$, then $\beta_{1}$ = 0. Moreover, if $\mu_{j}\neq 0$ for some $j=1,\ldots,d$, then the skewness index $\beta_{1}$ may diverge to $\infty$ or $-\infty$ as $p\rightarrow 0$. Therefore, for MGSN distribution Mardia’s multivariate skewness index varies from $-\infty$ to $\infty$.

If \displaystyle\mbox{\boldmathX}\sim MGSN{}_{d}(p,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}) and if we denote the mean vector and dispersion matrix of $X$, as {\mbox{\boldmath\mu}}_{\mbox{\boldmathX}} and {\mbox{\boldmath\Sigma}}_{\mbox{\boldmathX}}, respectively, then from (8), (9) and (10), we have the following relation:

[TABLE]

The following result provides the canonical correlation between {\mbox{\boldmathX}}_{1} and {\mbox{\boldmathX}}_{2}. It may be mentioned that canonical correlation is very useful in multivariate data analysis. In an experimental context suppose we take two sets of variables, then the canonical correlation can be used to see what is common among these two sets of variables; see for example Rao [12].

Theorem 2: Suppose {\mbox{\boldmathX}}\sim MGSN{}_{d}(p,{\mbox{\boldmath0}},{\mbox{\boldmath\Sigma}}). Further $X$ and $\Sigma$ are partitioned as in (7). Then for {\mbox{\boldmath\alpha}}\in\mathbb{R}^{h} and {\mbox{\boldmath\beta}}\in\mathbb{R}^{d-h} such that {\mbox{\boldmath\alpha}}^{T}{\mbox{\boldmath\Sigma}}_{11}{\mbox{\boldmath\alpha}}=1 and {\mbox{\boldmath\beta}}^{T}{\mbox{\boldmath\Sigma}}_{22}{\mbox{\boldmath\beta}}=1, the maximum corr({\mbox{\boldmath\alpha}}^{T}{\mbox{\boldmathX}}_{1},{\mbox{\boldmath\beta}}^{T}{\mbox{\boldmathX}}_{2}) = $\lambda_{1}$, where $\lambda_{1}$ is the maximum root of the $d$-degree polynomial equation

[TABLE]

Proof: From Theorem 1, we obtain

[TABLE]

Therefore, using (11), it follows that the problem is to find {\mbox{\boldmath\alpha}}\in\mathbb{R}^{h} and {\mbox{\boldmath\beta}}\in\mathbb{R}^{d-h} such that it maximizes

[TABLE]

subject to the restrictions {\mbox{\boldmath\alpha}}^{T}{\mbox{\boldmath\Sigma}}_{11}{\mbox{\boldmath\alpha}}=1 and {\mbox{\boldmath\beta}}^{T}{\mbox{\boldmath\Sigma}}_{22}{\mbox{\boldmath\beta}}=1. Now following the same steps as in the multivariate normal cases, Anderson [1], the result follows.

The following result provides the characteristic function of the Wishart type matrix based on MGSN random variables.

Theorem 3: Suppose $Z_{1},\ldots,Z_{n}$ are $n$ i.i.d. random variables, and $Z_{1}\sim$ MGSN{}_{d}(p,{\mbox{\boldmath0}},{\mbox{\boldmath\Sigma}}). Let us consider the Wishart type matrix

[TABLE]

If {\mbox{\boldmath\Theta}}=((\theta_{ij})) with $\theta_{ij}=\theta_{ji}$ is a $d\times d$ matrix, then the characteristic function of $(A_{11},\ldots,A_{pp}$, $2A_{12},2A_{13},\ldots,2A_{p-1,p})$ is

[TABLE]

Proof:

[TABLE]

Now we would like to compute \displaystyle E\left[e^{i\left({\mbox{\boldmathZ}}_{1}^{T}{\mbox{\boldmath\Theta}}{\mbox{\boldmathZ}}_{1}\right)}\right]. For a $d\times d$ real symmetric matrix $\Theta$, there is a real $d\times d$ matrix $B$ such that

[TABLE]

Here, $\displaystyle\hbox{diag}\{\delta_{1},\ldots,\delta_{d}\}$ means a $d\times d$ diagonal matrix with diagonal entries as $\delta_{1},\ldots,\delta_{d}$. If we make the transformation {\mbox{\boldmathZ}}_{1}={\mbox{\boldmathB}}{\mbox{\boldmathY}}, then using Theorem 1, {\mbox{\boldmathY}}\sim MGSN{}_{d}(p,{\mbox{\boldmath0}},{\mbox{\boldmathI}}). Using the definition MGSN distribution it follows that

[TABLE]

here $N\sim$ GE$(p)$, and {\mbox{\boldmathY}}_{m}’s are i.i.d. random vectors, and {\mbox{\boldmathY}}_{1}\sim N{}_{d}({\mbox{\boldmath0}},{\mbox{\boldmathI}}). We use the following notation

[TABLE]

Hence, $\displaystyle Y_{j}=\sum_{m=1}^{N}Y_{mj}$, for $j=1,\ldots,d$. Therefore,

[TABLE]

Theorem 4: If for any {\mbox{\boldmathc}}\neq{\mbox{\boldmath0}},{\mbox{\boldmathc}}\in\mathbb{R}^{d}, {\mbox{\boldmathc}}^{T}{\mbox{\boldmathX}}\sim GSN(\mu({\mbox{\boldmathc}}),\sigma({\mbox{\boldmathc}}),p) for a $d$ dimensional random vector $X$, then there exists a $d$-dimensional vector $\mu$ and a $d\times d$ symmetric matrix $\Sigma$ such that {\mbox{\boldmathX}}\sim MGSN{}_{d}(p,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}). Here -\infty<\mu({\mbox{\boldmathc}})<\infty,0<\sigma({\mbox{\boldmathc}})<\infty are functions of $c$.

Proof: If we denote the mean vector and dispersion matrix of the random vector $X$, as {\mbox{\boldmath\mu}}_{\mbox{\boldmathX}} and {\mbox{\boldmath\Sigma}}_{\mbox{\boldmathX}}, respectively, then we have E({\mbox{\boldmathc}}^{T}{\mbox{\boldmathX}})={\mbox{\boldmathc}}^{T}{\mbox{\boldmath\mu}}_{\small{\mbox{\boldmathX}}} and V({\mbox{\boldmathc}}^{T}{\mbox{\boldmathX}})={\mbox{\boldmathc}}^{T}{\mbox{\boldmath\Sigma}}_{\mbox{\boldmathX}}{\mbox{\boldmathc}}. Hence from (3) and (4), we have

[TABLE]

Therefore, from (1), using $t$ = 1, it follows that

[TABLE]

Let us define a $d$-dimensional vector $\mu$ and a $d\times d$ symmetric matrix $\Sigma$ as given below:

[TABLE]

Therefore,

[TABLE]

and (14) can be written as

[TABLE]

Hence the result follows.

Therefore, combining Theorems 1 and 2, we have the following characterization results for a $d$-variate MGSN distribution.

Theorem 5: If a $d$-dimensional random vector $X$ has a mean vector {\mu}_{\mbox{\boldmathX}} and a dispersion matrix {\mbox{\boldmath\Sigma}}_{\mbox{\boldmathX}}, then {\mbox{\boldmathX}}\sim MGSN{}_{d}(p,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}), here $\mu$ and $\Sigma$ are as defined in (15), if and only if for any {\mbox{\boldmathc}}\neq{\mbox{\boldmath0}},{\mbox{\boldmathc}}\in\mathbb{R}^{d}, {\mbox{\boldmathc}}^{T}{\mbox{\boldmathX}}\sim GSN(\mu({\mbox{\boldmathc}}),\sigma({\mbox{\boldmathc}}),p), where \mu({\mbox{\boldmathc}}) and \sigma({\mbox{\boldmathc}}) are as in (13).

Now we provide another characterization of the MGSN distribution.

Theorem 6: Suppose {\mbox{\boldmathX}}_{1},{\mbox{\boldmathX}}_{2},\ldots is a sequence of i.i.d. $d$-dimensional random vectors, and $M\sim$ GE$(\alpha)$, for $0<\alpha\leq 1$. Consider a new $d$-dimensional random vector

[TABLE]

Then {\mbox{\boldmathY}}\sim MGSN{}_{d}(\beta,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}) for $\beta\leq\alpha$, if and only if {\mbox{\boldmathX}}_{1} has a MGSN distribution.

Proof: If part. Suppose {\mbox{\boldmathX}}_{1}\sim MGSN{}_{d}(p,{\mbox{\boldmath\mu}},\mbox{\boldmath\Sigma}), then the MGF of $Y$ for {\mbox{\boldmatht}}\in\mathbb{R}^{d}, can be written as

[TABLE]

here $\beta=\alpha p\leq\alpha$.

Only if part. Suppose {\mbox{\boldmathY}}\sim MGSN{}_{d}(\beta,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}) for some $0<\beta\leq\alpha$, and the MGF of {\mbox{\boldmathX}}_{1} is M_{{\mbox{\boldmathX}}_{1}}({\mbox{\boldmatht}}). We have the following relation:

[TABLE]

From (16), we obtain

[TABLE]

for $\gamma=\beta/\alpha\leq 1$. Therefore, {\mbox{\boldmathX}}_{1}\sim MGSN{}_{d}(\gamma,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}).

Stochastic ordering plays a very important role in the distribution theory. It has been studied quite extensively in the statistical literature. For its importance and different applications, interested readers are referred to Shaked and Shantikumar [14]. Now we will discuss the multivariate total positivity of order two (MTP2) property, in the sense of Karlin and Rinott [8], of the joint PDF of MGSN distribution. We shall be using the following notation here. For any two real numbers $a$ and $b$, let $a\wedge b$ = min$\{a,b\}$, and $a\vee b$ = max$\{a,b\}$. For any vector {\mbox{\boldmathx}}=(x_{1},\ldots,x_{d})^{T} and {\mbox{\boldmathy}}=(y_{1},\ldots,y_{d})^{T}, let {\mbox{\boldmathx}}\vee{\mbox{\boldmathy}}=(x_{1}\vee y_{1},\ldots,x_{d}\vee y_{d})^{T} and {\mbox{\boldmathx}}\wedge{\mbox{\boldmathy}}=(x_{1}\wedge y_{1},\ldots,x_{d}\wedge y_{d})^{T}. Let us recall the definition of MTP2 property. A function $g:\mathbb{R}^{d}\rightarrow\mathbb{R}^{+}$ is said to have MTP2 property, in the sense of Karlin and Rinott [8], if \displaystyle g({\mbox{\boldmathx}})g({\mbox{\boldmathy}})\leq g({\mbox{\boldmathx}}\wedge{\mbox{\boldmathy}})g({\mbox{\boldmathx}}\vee{\mbox{\boldmathy}}), for all {\mbox{\boldmathx}},{\mbox{\boldmathy}}\in\mathbb{R}^{d}. We then have the following result for MGSN distribution.

Theorem 7: Let {\mbox{\boldmathX}}\sim MGSN{}_{d}(p,{\mbox{\boldmath0}},{\mbox{\boldmath\Sigma}}), and all the off-diagonal elements of {\mbox{\boldmath\Sigma}}^{-1} are less than or equal to zero, then the PDF of $X$ has MTP2 property.

Proof: To prove that the PDF of $X$ has MTP2 property, it is enough to show that

[TABLE]

for any {\mbox{\boldmathx}}=(x_{1},\ldots,x_{d})^{T}\in\mathbb{R}^{d} and {\mbox{\boldmathy}}=(y_{1},\ldots,y_{d})^{T}\in\mathbb{R}^{d}. If the elements of {\mbox{\boldmath\Sigma}}^{-1} are denoted by $((\sigma^{kj}))$, for $k,j=1,\ldots,d$, then proving (17) is equivalent to showing

[TABLE]

For all $k,j=1,\ldots,d$,

[TABLE]

which can be easily shown by taking any ordering of $x_{k},x_{j},y_{k},y_{j}$. Now the result follows since $\sigma^{jk}\leq 0$.

The following two decompositions of a MGSN distribution are possible. We use the following notations. The distribution of a negative binomial random variable with parameters $r$ and $p$, where $r$ is a non-negative integer and $0<p<1$, will be denoted by NB$(r,p)$. If $T\sim$ NB$(r,p)$, then the MGF of $T$ is

[TABLE]

A discrete random variable $Z$ is said to have a logarithmic distribution with parameter $p$, for $0<p<1$, if the PMF of $Z$ is

[TABLE]

and it will be denoted by LD$(p)$. Now we provide two decompositions of MGSN distribution.

Decomposition 1: Suppose {\mbox{\boldmathX}}\sim MGSN{}_{d}(p,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}). Further, for any positive integer $n$ and for any $1\leq k\leq n$, suppose

[TABLE]

where $T\sim$ NB$(r,p)$ and $r=1/n$, {\mbox{\boldmathY}}_{j}’s are i.i.d. N{}_{d}(r{\mbox{\boldmath\mu}},r{\mbox{\boldmath\Sigma}}), $T$ and {\mbox{\boldmathY}}_{j}’s are independently distributed, then

[TABLE]

Proof: The MGF of {\mbox{\boldmathZ}}_{kn} can be written as

[TABLE]

It implies that MGSN law is infinitely divisible. The following decomposition is also possible.

Decomposition 2: Suppose $Q$ is a Poisson random variable with parameter $\lambda$, and $\{Z_{i},i=1,2,\ldots\}$ is a sequence of i.i.d. random variables having logarithmic distribution with the following probability mass function for $\lambda=-\ln p$;

[TABLE]

and all the random variables are independently distributed. If {\mbox{\boldmathX}}\sim MGSN{}_{d}(p,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}), then the following decomposition is possible

[TABLE]

here \{{\mbox{\boldmathY}}_{i}|Z_{i}=k\}\sim N{}_{d}(k{\mbox{\boldmath\mu}},k{\mbox{\boldmath\Sigma}}) for $i=1,2,\ldots$, and they are independently distributed, {\mbox{\boldmathY}}\sim N{}_{d}({\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}), and it is independent of $Q$, and ({\mbox{\boldmathY}}_{i},Z_{i}) for all $i=1,2,\ldots.$

Proof: First note that the probability generating function of $Q$ and $Z_{1}$ are as follows:

[TABLE]

The MGF of {\mbox{\boldmathY}}_{i} for {\mbox{\boldmatht}}\in\mathbb{R}^{d}, such that \displaystyle(1-p)e^{{\mbox{\boldmath\mu}}^{T}{\mbox{\boldmatht}}+\frac{1}{2}{\mbox{\boldmatht}}^{T}{\mbox{\boldmath\Sigma}}{\mbox{\boldmatht}}}<1, can be obtained as

[TABLE]

Therefore, the MGF of the right hand side of (18) can be written as

[TABLE]

The following results will be useful for further development. Let us consider the random vector ({\mbox{\boldmathX}},N), where $X$ and $N$ are same as defined in (5). The joint PDF of ({\mbox{\boldmathX}},N) can be written as

[TABLE]

for {\mbox{\boldmathx}}\in\mathbb{R}^{d} and for any positive integer $n$. Therefore, the conditional probability mass function of $N$ given {\mbox{\boldmathX}}={\mbox{\boldmathx}} becomes

[TABLE]

Therefore,

[TABLE]

and

[TABLE]

4 Statistical Inference

4.1 Estimation

In this section we discuss the maximum likelihood estimators (MLEs) of the unknown parameters, when $0<p<1$. When $p$ = 1, the MLEs of $\mu$ and $\Sigma$ can be easily obtained as the sample mean and the sample variance covariance matrix, respectively. Suppose {\cal D}=\{{\mbox{\boldmathx}}_{1},\ldots,{\mbox{\boldmathx}}_{n}\} is a random sample of size $n$ from MGSN{}_{d}(p,{\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}), then the log-likelihood function becomes

[TABLE]

The maximum likelihood estimators (MLEs) of the unknown parameters can be obtained by maximizing (21) with respect to the unknown parameters. It involves solving a $(d+1+d(d+1)/2)$ dimensional optimization problem. Therefore, for large $d$, it is a challenging issue.

To avoid that problem, first it is assumed that $p$ is known. For a known $p$, we estimate the MLEs of $\mu$ and $\Sigma$ by using EM algorithm, say \widehat{\mbox{\boldmath\mu}}(p) and \widehat{\mbox{\boldmath\Sigma}}(p), respectively. We maximize l(p,\widehat{\mbox{\boldmath\mu}}(p),\widehat{\mbox{\boldmath\Sigma}}(p)) to compute $\widehat{p}$, the MLE of $p$. Finally we obtain the MLE of $\mu$ and $\Sigma$ as \widehat{\mbox{\boldmath\mu}}=\widehat{\mbox{\boldmath\mu}}(\widehat{p}) and \widehat{\mbox{\boldmath\Sigma}}=\widehat{\mbox{\boldmath\Sigma}}(\widehat{p}), respectively. Now we will show how to compute \widehat{\mbox{\boldmath\mu}}(p) and \widehat{\mbox{\boldmath\Sigma}}(p), for a given $p$ using EM algorithm. We treat the problem as a missing value problem, and the main idea is as follows.

It is assumed that $p$ is known. Suppose we have the complete observations of the form \{({\mbox{\boldmathx}}_{1},m_{1}),\ldots,({\mbox{\boldmathx}}_{n},m_{n})\} from ({\mbox{\boldmathX}},N). Then the log-likelihood function based on the complete observation becomes (without the additive constant)

[TABLE]

Therefore, if we define the MLEs of $\mu$ and $\Sigma$ based on the complete observations as \widehat{\mbox{\boldmath\mu}}_{c}(p) and $\widehat{\Sigma}_{c}(p)$, respectively, then for $\displaystyle K=\sum_{i=1}^{n}m_{i}$,

[TABLE]

and

[TABLE]

Note that \displaystyle\widehat{\mbox{\boldmath\mu}}_{c}(p) is obtained by taking derivative of \displaystyle\sum_{i=1}^{n}\frac{1}{m_{i}}({\mbox{\boldmathx}}_{i}-m_{i}{\mbox{\boldmath\mu}})^{T}{\mbox{\boldmath\Sigma}}^{-1}({\mbox{\boldmathx}}_{i}-m_{i}{\mbox{\boldmath\mu}}) with respect to $\mu$, and equate it to zero. Similarly, \displaystyle\widehat{\mbox{\boldmath\Sigma}}_{c}(p) is obtained by using Lemma 3.2.2 of Anderson [1].

Now we are ready to provide the EM algorithm for a given $p$. The EM algorithm consists of maximizing the conditional expectation of the complete log-likelihood function, based on the observed data and the current value of {\mbox{\boldmath\theta}}=({\mbox{\boldmath\mu}},{\mbox{\boldmath\Sigma}}), say \widetilde{\mbox{\boldmath\theta}}, in an iterative two-step algorithm process, see for example Dempster et al. [6]. The E-step is to compute the conditional expectation denoted by Q({\mbox{\boldmath\theta}}|\widetilde{\mbox{\boldmath\theta}}), and the M-step is maximizing Q({\mbox{\boldmath\theta}}|\widetilde{\mbox{\boldmath\theta}}), with respect to $\theta$. We use the following notations:

[TABLE]

where $a_{i}$ and $b_{i}$ are obtained using (19) and (20), respectively.

E-Step: It consists of calculating Q({\mbox{\boldmath\theta}}|\widetilde{\mbox{\boldmath\theta}}), \widetilde{\mbox{\boldmath\theta}} being the current parameter value.

[TABLE]

M-Step: It involves maximizing Q({\mbox{\boldmath\theta}}|\widetilde{\mbox{\boldmath\theta}}) with respect to $\theta$, to obtain \overline{\mbox{\boldmath\theta}}, where

[TABLE]

Here, \displaystyle\hbox{arg max}_{\mbox{\boldmath\theta}}Q({\mbox{\boldmath\theta}}|\widetilde{\mbox{\boldmath\theta}}) means the value of $\theta$ for which the function \displaystyle Q({\mbox{\boldmath\theta}}|\widetilde{\mbox{\boldmath\theta}}) takes the maximum value. From (22) and (23), we obtain

[TABLE]

and

[TABLE]

We propose the following algorithm to compute the MLEs of $\mu$ and $\Sigma$ for a known $p$.

Algorithm 2:

• •

Step 1: Choose an initial guess of $\theta$, say {\mbox{\boldmath\theta}}^{(0)}.

• •

Step 2: Obtain

[TABLE]

• •

Step 3: Continue the process until convergence takes place.

Once for a given $p$, the MLEs of $\mu$ and $\Sigma$ are obtained, say \widehat{\mbox{\boldmath\mu}}(p) and \widehat{\mbox{\boldmath\Sigma}}(p), respectively, then the MLE of $p$ can be obtained by maximizing the profile log-likelihood function of $p$, i.e. \displaystyle l(p,\widehat{\mbox{\boldmath\mu}}(p),\widehat{\mbox{\boldmath\Sigma}}(p)), with respect to $p$. If it is denoted by $\widehat{p}$, then the MLEs of $\mu$ and $\Sigma$ become \widehat{\mbox{\boldmath\mu}}=\widehat{\mbox{\boldmath\mu}}(\widehat{p}) and \widehat{\mbox{\boldmath\Sigma}}=\widehat{\mbox{\boldmath\Sigma}}(\widehat{p}), respectively. The details will be explained in Section 5. We have used the sample mean vector and the sample variance covariance matrix as the initial guess of $\mu$ and $\Sigma$, respectively, of the proposed EM algorithm for all $p$.

4.2 Testing of Hypotheses

In this section we discuss three different testing of hypotheses problems which can be useful in practice. We propose to use the likelihood ratio test (LRT) in all the cases, and we indicate the asymptotic distribution of the LRT tests under the null hypothesis in each case. With the abuse of notations, in each case if $\delta$ is any unknown parameter, the MLE of $\delta$ under the null hypothesis will be denoted by $\widehat{\delta}_{H}$.

Test 1:

[TABLE]

The above testing problem (25) is important in practice as it tests the normality of the distribution. In this \widehat{\mbox{\boldmath\mu}}_{H} and \widehat{\mbox{\boldmath\Sigma}}_{H}, respectively, can be obtained as the sample mean and the sample variance covariance matrix. Since in this case $p$ is in the boundary, the standard results do not work. But using result 3 of Self and Liang [13] it follows that under the null hypothesis

[TABLE]

Test 2:

[TABLE]

The above testing problem (26) is important as it tests the symmetry of the distribution. In this case under the null hypothesis the MLEs of $p$ and $\Sigma$ can be obtained as follows. For a given $p$, the MLE of $\Sigma$ can be obtained using the EM algorithm as before, and then the MLE of $p$ can be obtained by maximizing the profile likelihood function. In this case the ’E-step’ and ’M-Step’ can be obtained from (4.1) and (24), respectively, by replacing {\mbox{\boldmath\mu}}={\mbox{\boldmath0}}. Under $H_{0}$, then

[TABLE]

Test 3:

[TABLE]

The above testing problem (27) is important as it tests the uncorrelatedness of the components. In this case the diagonal elements of the matrix $\Sigma$ will be denoted by $\sigma_{1}^{2},\ldots,\sigma_{d}^{2}$, i.e. \displaystyle{\mbox{\boldmath\Sigma}}=\hbox{diag}\left\{\sigma_{1}^{2},\ldots,\sigma_{d}^{2}\right\}. Now we will mention how to compute the MLEs of the unknown parameters $p$, $\mu$ and $\sigma_{1}^{2},\ldots,\sigma_{d}^{2}$, under the null hypothesis. In this case also as before for a given $p$, we use EM algorithm to compute the MLEs of $\mu$ and $\sigma_{1}^{2},\ldots,\sigma_{d}^{2}$, and finally the MLE of $p$ can be obtained by maximizing the profile likelihood function. Now we will describe how to compute the MLEs of $\mu$ and $\sigma_{1}^{2},\ldots,\sigma_{d}^{2}$, for a given $p$, by using the EM algorithm. We use the following notation for further development. The matrix \displaystyle{\mbox{\boldmath\Delta}}_{k} is a $d\times d$ matrix with all the entries 0, except the $(k,k)$-th element which is 1. Now under $H_{0}$, the ‘E-Step’ of the EM algorithm can be written as follows:

[TABLE]

The ‘M-Step’ involves maximizing (28) with respect to $\mu$, $\sigma_{1}^{2},\ldots,\sigma_{d}^{2}$ to obtain updated $\mu$, $\sigma_{1}^{2},\ldots,\sigma_{d}^{2}$, say \overline{\mbox{\boldmath\mu}}, \overline{\mbox{\boldmath\sigma}^{2}_{1}},\ldots,\overline{\mbox{\boldmath\sigma}^{2}_{d}}, respectively. From (28), we obtain

[TABLE]

and

[TABLE]

Here for a square matrix $A$, {\mbox{\boldmathA}}_{k} denotes the $k$-th diagonal element of the matrix $A$. Under the null hypothesis

[TABLE]

5 Simulations and Data Analysis

In this section we perform some Monte Carlo simulations to show how the proposed EM algorithm performs and we perform the analyses of two data sets analysis to show how the proposed model and the methods can be used in practice.

5.1 Simulation Results

For simulation purposes we have used the following sample size and the parameter values;

[TABLE]

Now to show the effectiveness of the EM algorithm we have considered both the cases namely when (a) $p$ is known and (b) $p$ is unknown. We have generated samples from the above configuration and computed the MLEs of $\mu$ and $\Sigma$ using EM algorithm. In all the cases we have used the sample mean and the sample variance covariance matrix as the initial guesses of the EM algorithm. We replicate the process 1000 times and report the average estimates and the associated mean squared errors (MSEs). For known $p$, the results are reported in Tables 1 and 3 and for unknown $p$, the results are reported in Tables 2 and 4. In each box of a table, the first figure, second figure and the third figure represent the true value, the average estimate and the corresponding MSE, respectively.

It is clear that the performances of the proposed EM algorithm are quite satisfactory. It is observed that the sample mean and the sample variance covariance matrix can be used as good initial guesses of the EM algorithm. In all the cases considered it is observed that the EM algorithm converges within 30 iterations, hence it can be used in practice quite conveniently. Further, it is observed that the profile likelihood method is also quite effective in estimating $p$, when it is unknown.

In this section we present the analysis of two data sets namely (i) one simulated data set and (ii) one real data set mainly to illustrate how the proposed EM algorithm and the other testing procedures can be used in practice.

5.2 Simulated Data Set

We have generated a data set using the Algorithm 1 as suggested in Section 2, with the following specification:

[TABLE]

It is available in http://home.iitk.ac.in/$\sim$kundu/fort.76. We present some basic statistics of the data set. The sample mean vector, and the sample variance covariance matrix are as follows:

[TABLE]

We start the EM algorithm for each $p$ with the above initial guesses. The profile log-likelihood function is plotted in Figure 2.

Finally, the MLEs of the unknown parameters are obtained as follows:

[TABLE]

and the associated log-likelihood value is -741.347. It may be mentioned for each $p$, the EM algorithm is continued for 20 iterations, and the log-likelihood value (21) is calculated based on the first 50 terms of the infinite series. The program is written in FORTRAN-77, and it is available in http://home.iitk.ac.in/$\sim$kundu/mv-geo-sn-em-punknown-data.for.

For illustrative purposes, we would like to perform the test:

[TABLE]

Under $H_{0}$, the MLEs of $\mu$ and $\Sigma$ become \bar{\mbox{\boldmathx}} and $S$, respectively, as given in (29), and the associated log-likelihood value is -887.852. Therefore, the value of the test statistic $T_{1}$ = 293.01, and the associated $p$ value is less than 0.00001. Hence, we reject $H_{0}$. Next we consider the following testing problem

[TABLE]

In this case under $H_{0}$, the MLEs of $p$ and $\Sigma$ are as follows

[TABLE]

and the associated log-likelihood value is -917.674. In this case the value of the test statistics $T_{2}$ = 352.654. Since the associated $p$ value is less than 0.00001, we reject the null hypothesis. Finally we consider the testing problem:

[TABLE]

In this case under $H_{0}$, the MLEs of the unknown parameters are as follows:

[TABLE]

and

[TABLE]

The associated log-likelihood value is -1036.80. The value of $T_{3}$ = 590.91. In this case also we reject $H_{0}$, as the associated $p$ values is less than 0.00001.

5.3 Stiffness Data Set

In this section we present the analysis of a real data set to show how the proposed model and the methodologies work in practice. The data set represents the four different measurements of stiffness, $x_{1},x_{2},x_{3},x_{4}$ of ‘Shock’ and ‘Vibration’ of each of 30 boards. The first measurement (Shock) involves sending a shock wave down the board and the second measurement (Vibration) is determined while vibrating the board. The last two measurements are obtained from static tests. The data set is available in Johnson and Wichern [7]. For easy reference it is presented in Table 5. Since all the entries of the data set are non-negative, if we want to the fit the multivariate skew normal distribution to this data set, the MLEs of the unknown parameters may not exist. In fact we have tried to fit univariate skew-normal distribution to $x_{1}$ and it is observed that the likelihood function is an increasing function of the ‘tilt’ parameter for fixed location and scale parameters. Therefore, the MLEs do not exist in this case. It is expected the same phenomenon even for skew-$t$ distribution for large values of the degrees of freedom.

Before progressing further we have divided all the measurements by 100, and it is not going to make any difference in the inferential procedure. The sample mean vector and the sample variance covariance matrix of the transformed data are

[TABLE]

Based on the EM algorithm and using the profile likelihood method, we obtain the MLEs of the unknown parameters as follows:

[TABLE]

and the associated log-likelihood value is -271.969. Now to check whether the proposed MGSN distribution provides a better fit than the multivariate normal distribution or not, we perform the following test:

[TABLE]

Under $H_{0}$, the MLEs of $\mu$ and $\Sigma$ are provided in (30), and the associated log-likelihood value is -277.761. Therefore, the value of the test statistic $T_{1}$ = 11.584, and the associated $p$ value is 0.0000025. Hence, we reject the null hypothesis, and it indicates that the proposed MGSN distribution provides a better fit than the multivariate normal distribution to the given stiffness data set. AIC also prefers MGSN distribution than the multivariate normal distribution for this data set.

6 Conclusion

In this paper we have discusses different properties of the MGSN distribution in details. Different characterization results and dependence properties have been established. The $d$-dimensional MGSN distribution has $d+1+d(d+1)/2$ unknown parameters. We have proposed to use EM algorithm and the profile likelihood method to compute the MLEs of the unknown parameters, and it is observed that the proposed algorithm can be implemented very easily. We have discussed some testing of hypothesis problems also. Two data sets have been analyzed to show the effectiveness of the proposed methods, and it is observed that for the real ’stiffness’ data set MGSN provides a better fit than the multivariate normal distribution. Hence, this model can be used as an alternative to Azzalini’s multivariate skew normal distribution.

Acknowledgements

The author would like to thank Prof. Dimitris Karlis from the Athens University of Economics and Business, Greece, for his helpful comments. The author is really thankful to two reviewers and the associate editor for their careful reading and providing constructive comments.