# Moments and random number generation for the truncated elliptical family of distributions

This paper proposes an algorithm to generate random numbers from any member of the truncated multivariate elliptical family of distributions with a strictly decreasing density generating function. Based on Neal (2003) and Ho et al. (2012), we construct an efficient sampling method by means of a slice sampling algorithm with Gibbs sampler steps. We also provide a faster approach to approximate the first and the second moment for the truncated multivariate elliptical distributions where Monte Carlo integration is used for the truncated partition, and explicit expressions for the non-truncated part (Galarza et al., 2020). Examples and an application to environmental spatial data illustrate its usefulness. Methods are available for free in the new R library elliptical.

## Authors

• 2 publications
• 6 publications
• 8 publications
09/28/2020

### On moments of folded and doubly truncated multivariate extended skew-normal distributions

This paper develops recurrence relations for integrals that relate the d...
04/04/2022

### Scalable random number generation for truncated log-concave distributions

Inverse transform sampling is an exceptionally general method to generat...
04/20/2022

### Integral, mean and covariance of the simplex-truncated multivariate normal distribution

Compositional data, which is data consisting of fractions or probabiliti...
03/02/2022

### Doubly truncated moment risk measures for elliptical distributions

In this paper, we define doubly truncated moment (DTM), doubly truncated...
03/29/2022

### Direct Sampling with a Step Function

The direct sampling method proposed by Walker et al. (JCGS 2011) can gen...
08/09/2021

### Identification in Bayesian Estimation of the Skewness Matrix in a Multivariate Skew-Elliptical Distribution

Harvey et al. (2010) extended the Bayesian estimation method by Sahu et ...
10/20/2017

### Nonparametric estimation of multivariate distribution function for truncated and censored lifetime data

In this article we consider a number of models for the statistical data ...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

The use of truncated distributions arises in a wide variety of statistical models as survival analysis, censored data models, Bayesian models with truncated parameters space, and abound in such fields as agronomy, biology, environmental monitoring, medicine, and economics. Algorithms like Expectation-Maximization (EM)

(dempster1977maximum) are employed frequently in multivariate censored data analysis under a likelihood-based perspective for its facility to deal with missing and partially observed data. This algorithm requires the computation of conditional truncated moments, commonly the first two moments. For example, matos2013likelihood and matos2016censoredestimated the parameters of a censored mixed-effects model for irregularly repeated measures via the EM algorithm, which needed to compute the first two moments of a truncated multivariate (TMVT) and a truncated multivariate normal (TMVN) distributions, respectively.

In this context, there are a few libraries in R (Rproj) which provide truncated multivariate moments. For instance, the package tmvtnorm (wilhelm2015package)

computes the mean and the variance of the TMVN distribution by deriving its moment generating function, while the

MomTrunc library (galarza2021package) uses a recursive approach method proposed by kan2017moments to compute arbitrary higher-order moments. On the other hand, for the TMVT distribution, the packages TTmoment (ho2015r) and MomTrunc

compute its two first moments. Moreover, the first library only handles integer degrees of freedom greater than 4, while the latter can compute even high-order moments for any degrees of freedom

(galarza2021moments).

Variations of the EM algorithm such as Stochastic Approximation EM (SAEM) (delyon1999convergence) and Monte Carlo EM (MCEM) (wei1990monte) replace the conditional expectations by an approximation that requires to draw independent random observations from a truncated distribution. For instance, lachos2017influence estimated the parameters of a linear spatial model for censored data using the SAEM algorithm, which needed to generate random samples from the TMVN distribution to perform the stochastic approximation step. More recently, also using the SAEM algorithm, lachos2019flexible

proposed a robust multivariate linear mixed model for multiple censored responses based on the scale mixtures of normal (SMN) distributions. Moreover, generating random numbers from truncated distributions is also required in Bayesian models,

gelfand1992bayesian showed how to perform Bayesian analysis for constrained parameters or truncated data problems by using Gibbs sampling.

There are several methods to generate random samples from a truncated distribution in the literature, and the common one is the rejection sampling technique. This method draws samples from the non-truncated distribution and retains only the samples inside the support region. However, the procedure may be inefficient, especially when the truncation interval is too small or it is located at a less probable area of the probability density function (pdf).

neal2003slice

proposed proposed the Slice sampling method, a procedure that turns sampling from a truncated density into sampling repeatedly from uniform distributions instead. This algorithm is easy to code, fast and does not reject samples, making it more efficient than the conventional rejection method.

To the best of our knowledge, there are no proposals in the literature to generate samples from other multivariate truncated distributions in the elliptical class other than the TMVN and TMVT distributions (available in the tmvtnorm and TTmoment packages). Hence, motivated by the slice sampling algorithm, we propose a general method to obtain samples from any truncated multivariate elliptical distribution with strictly decreasing density generating function (dgf). Using conditional expectation properties, we also propose an efficient algorithm to approximate the moments of the most common distribution of this class: the truncated multivariate normal, Student-

, slash, contaminated normal, and Pearson VII distributions. This method requires less running time when compared with the existing ones, since it deals with the truncated and non-truncated part of the vector separately. Our proposal can be reached through the

R package relliptical

. Finally, it is worth mentioning that moments of truncated elliptical distributions can be used to compute truncated moments for the selection elliptical family of distributions, a wide family which includes complex multivariate asymmetric versions of the elliptical distributions as the extended skew-normal, the unified skew-

distributions, among others. Therefore, our proposal opens the doors for the calculation of truncated moments of complex elliptical asymmetric distributions, which are of particular interest for the development of robust censored models with asymmetry, heavy tails and missingness (see for instance galarza2021skew; de2021finite).

The paper is organized as follows. Section 2 shows some results related to the elliptical and truncated elliptical family of distributions and a brief description of the slice sampling algorithm. Section 3 is devoted to the formulation of the sampling algorithm for the truncated elliptical distributions, whereas Section 4 focuses on our proposed method to approximate the first and the second moment. For the last two sections, we present a brief introduction to its respective R function. A simulation study that compares the mean and covariance matrix for the TMVT distribution estimated through different methods in R is presented as well. Section 5 displays an application on censored Gaussian spatial models throughout the analysis of the Missouri dioxin contamination dataset. Finally, Section 6 concludes with a discussion.

## 2 Preliminaries

### 2.1 Elliptical Family of Distributions

As defined in muirhead2009aspects and fang2018symmetric, a random vector is said to follow an elliptical distribution with location parameter , positive-definite scale matrix , and density generating function , if its pdf is given by

 fX(x)=cp|Σ|−1/2g((x−μ)⊤Σ−1(x−μ)),x∈Rp, (1)

where is a non-negative Lebesgue measurable function on such that and denotes the determinant of matrix . Moreover,

 cp=Γ(p/2)πp/2(∫∞0tp/2−1g(t)dt)−1

is the normalizing constant, with representing the complete gamma function. We will use the notation .

Members of the elliptical family of distributions are characterized by their density generating function . Some examples of the elliptical family of distributions are:

• [leftmargin=.4cm]

• The multivariate normal distribution, , with mean and variance-covariance matrix , arises when the dgf takes the form .

• The multivariate Student-t distribution, , where is the location parameter, is the scale matrix, and is called the degrees of freedom, is obtained when .

• The multivariate power exponential,

, with kurtosis parameter

. In this case,

. A particular case of the power exponential distribution is the normal distribution, which arises when

.

• The multivariate slash,

, we get a random variable with multivariate slash distribution when

.

• The multivariate Pearson VII distribution, , with parameters , , , and is obtained when .

For more distributions belonging to this family, please see fang2018symmetric.

### 2.2 Truncated Elliptical Family of Distributions

Let be a measurable set. We say that a random vector has truncated elliptical distribution with support , location parameter , scale parameter and dgf , if its pdf is given by

 fY(y)=g((y−μ)⊤Σ−1(y−μ))∫Ag((y−μ)⊤Σ−1(y−μ))dy=fX(y)% Pr(X∈A),y∈A, (2)

where . We use the notation . Notice that the pdf of Y is written as the ratio between the pdf of and , so the pdf of Y exists if the pdf of X does, which occurs if is a positive-definite (see, moran2019new). The variable Y is also said to be an elliptical distribution truncated on , being represented by .

As in the elliptical family of distributions, the dgf determines any distribution within the truncated elliptical class of distributions, for example, if , then Y has TMVT distribution. We will denote the different members of the truncated elliptical family defined in the subsection before as for the TMVN distribution, for the TMVT distribution, for the truncated multivariate power exponential, for the truncated multivariate slash distribution, and for the truncated multivariate Pearson VII distribution.

### 2.3 Slice Sampling Algorithm

Introduced by neal2003slice

, the slice sampling algorithm is a Markov Chain Monte Carlo (MCMC) method for drawing random samples from a given distribution. The algorithm’s idea is to sample uniformly from the

-dimensional region under the graph of , a non-negative function proportional to the pdf of X. Hence, let be an auxiliary variable such that the joint pdf of X and is uniform over the region , i.e., , with being the indicator function. Therefore, we can obtain samples from the distribution of X by sampling jointly and then ignoring values.

Note that generating independent random points uniformly distributed on may not be easy. To overcome this problem, neal2003slice defined a Markov Chain that converges to an uniform distribution, in the same manner than the Gibbs sampling or Metropolis-Hastings algorithms. Then, considering Gibbs sampler steps, the slice sampling algorithm at iteration works as follows: given the current value of sample from , then draw from the conditional distribution of X given , which is uniform over the region , i.e., , for all , where is the desired sample size.

Figure 1 shows the steps of the slice sampling algorithm for being a univariate random variable. Given an initial value , we draw uniformly over the interval and then we sample from the conditional distribution of , i.e., uniformly over the interval . These two steps are repeated times.

## 3 Sampling from the Truncated Elliptical Family of Distributions

Next, we describe the proposed slice sampling algorithm with Gibbs sampler steps to generate samples from a multivariate elliptical distribution with strictly decreasing dgf. Without loss of generality, we first consider a -variate truncated elliptical distribution with zero location parameter, positive-definite scale matrix , dgf , and truncation region , , in other words, we will consider . Here is a correlation matrix, such that the scale matrix can be written as , where . The pdf of X is given by

 fX(x)∝g(x⊤R−1x)I(x∈A), (3)

Now, in order to sample uniformly from the -dimensional region under the plot of , we introduce an auxiliary variable , such that the joint pdf of X and is

 fX,Y(x,y)∝I(0

It is enough to calculate the conditional distributions of and in order to established our slice sampling algorithm with Gibbs steps to generate independent random observations from the pdf in (4). These are given by:

 fY|X(y|x) ∝ I(0

Note that sampling from the distribution of is straightforward, but sampling from is not trivial. Thus, we use the idea of ho2012some, that consists in sampling each element of X given the remaining elements, i.e., sampling given and , for all . Hence, the following steps are performed to draw a random number from the distribution of .

1. [leftmargin=0.6cm]

2. Let . Since is a strictly decreasing function, it follows that is equivalent to .

3. Write , where is the th element of the inverse of R, and .

4. Combining items 1 and 2, we obtain that , where

5. Because , thereby .

Therefore, the steps to draw samples from a -variate truncated elliptical distribution are summarized in Algorithm 1. As seen, only univariate uniform simulations are involved in the algorithm which are fast to compute. Note also that the assumption that the dgf is strictly decreasing has been used in step 1. A general case can be easily considered by studying the extrema points of .

Moreover, members of the truncated elliptical family of distributions are closed under affine transformations (fang2018symmetric). Hence drawing samples from may be readily done by sampling first from , and then recovery Y by the following transformation , such that , , and .

### 3.1 R function and Examples

Algorithm 1 and the transformation described previously were implemented in the R package relliptical. Its main function for random number generation is called rtelliptical, whose signature is the following.

rtelliptical(n=1e4, mu=rep(0,length(lower)), Sigma=diag(length(lower)), lower,
upper=rep(Inf,length(lower)), dist="Normal", nu=NULL, expr=NULL,
gFun=NULL, ginvFun=NULL, burn.in=0, thinning=1)

In this function, is the number of observations to be sampled, nu is the additional parameter or vector of parameters depending on the distribution of X, mu is the location parameter, Sigma is the positive-definite scale matrix, and lower and upper are the lower and upper truncation points, respectively. The truncated normal, Student-, power exponential, Pearson VII, slash, and contaminated normal distributions can be specified through the argument dist.

The following examples illustrate the function rtelliptical, for drawing samples from truncated bivariate distributions with location parameter , scale matrix elements , and , and truncation region , with and . The distributions considered are the predefined ones in the package.

• [leftmargin=.4cm]

• Truncated normal

rtelliptical(n=1e4, mu=c(0,0), Sigma=matrix(c(1,0.7,0.7,1),2,2), lower=c(-2,-2),
upper=c(3,2), dist="Normal")
• Truncated Student- with degrees of freedom

rtelliptical(n=1e4, mu=c(0,0), Sigma=matrix(c(1,0.7,0.7,1),2,2), lower=c(-2,-2),
upper=c(3,2), dist="t", nu=3)
• Truncated power exponential with kurtosis

rtelliptical(n=1e4, mu=c(0,0), Sigma=matrix(c(1,0.7,0.7,1),2,2), lower=c(-2,-2),
upper=c(3,2), dist="PE", nu=2)
• Truncated Pearson VII with parameters and

rtelliptical(n=1e4, mu=c(0,0), Sigma=matrix(c(1,0.7,0.7,1),2,2), lower=c(-2,-2),
upper=c(3,2), dist="PVII", nu=c(2.50, 3.0))
• Truncated slash with 3/2 degrees of freedom

rtelliptical(n=1e4, mu=c(0,0), Sigma=matrix(c(1,0.7,0.7,1),2,2), lower=c(-2,-2),
upper=c(3,2), dist="Slash", nu=1.50)
• Truncated contaminated normal with and

rtelliptical(n=1e4, mu=c(0,0), Sigma=matrix(c(1,0.7,0.7,1),2,2), lower=c(-2,-2),
upper=c(3,2), dist="CN", nu=c(0.70, 0.20))

Note that, no additional arguments are passed for the TMVN distribution. In the opposite way, for the truncated contaminated normal and Pearson VII distributions, nu is a vector of length two, and for the remaining distributions, this parameter is a non-negative scalar. An important remark is that exists closed form expressions to compute for the normal, Student-, power exponential, and Pearson VII distributions, however, the contaminated normal and slash distributions require numerical methods for this purpose. This value is calculated as the root of the function , through the Newton-Raphson algorithm for the contaminated normal, and using Brent’s method (brent2013algorithms)

, for the slash distribution, a mixture of linear interpolation, inverse quadratic interpolation, and the bisection method.

This function also allows generating random numbers from other truncated elliptical distributions not specified in the dist argument, by supplying the dgf through arguments either expr or gFun. The easiest way is to provide the dgf expression to argument expr as a character. The notation used in expr needs to be understood by package Ryacas0 (andersen2020ryacas), and the R environment. For instance, for the dgf , the user must provide expr = "exp(1)^(-t)". For this case, when a character expression is provided to expr, the algorithm tries to compute a closed-form expression for the inverse function of , however, this is not always possible (a warning message is returned). On the other hand, if it is no possible to pass an expression to expr, due to the complexity of the expression, the user may provide a custom R function to the gFun argument. By default, its inverse function is approximated numerically, however, the user may also provide its inverse to the ginvFun argument to gain some computational time. When gFun is provided, arguments dist and expr are ignored.

For example, to generate samples from the bivariate truncated logistic distribution with same parameters as before, and which has dgf , we can run the following code.

rtelliptical(n=1e4, mu=c(0,0), Sigma=matrix(c(1,0.7,0.7,1),2,2), lower=c(-2,-2),
upper=c(3,2), expr="exp(1)^(-t)/(1+exp(1)^(-t))^2")

Another distribution that belongs to the elliptical family is the Kotz-type distribution with parameters , and , whose dgf is (fang2018symmetric). For this distribution, is not strictly decreasing, however, for , it holds. Hence, our proposal works for , , and . For this type of more complex dgf, it is advisable to pass it through the gFun argument as an R function (with other parameters as fixed values). In the following example, we draw samples from a bivariate Kotz-type distribution with settings as before, and extra parameters , and .

rtelliptical(n=1e4, mu=c(0,0), Sigma=matrix(c(1,0.7,0.7,1),2,2), lower=c(-2,-2),
upper=c(3,2), gFun=function(t){t^(-1/2)*exp(-2*t^(1/4))})

Figure 2 shows the scatterplot and marginal histograms for the observations sampled from each of the truncated bivariate distributions referred above.

As mentioned by robert2010introducing and ho2012some, the slice sampling algorithm with Gibbs steps generates random samples conditioned on previous values, resulting in a sequence of correlated samples. Thus, it is essential to analyze the dependence effect of the proposed algorithm. Figure 6 in Section A.1 displays the autocorrelation plots for each one of the distributions, where we notice that the autocorrelation drops quickly and becomes negligibly small when lags become large, evidencing well mixing and quickly converging for these examples. If necessary, initial observations can be discarded by means of the burn.in argument. Finally, autocorrelation can be decimated by setting the thinning argument. Thinning consists in picking separated points from the sample, at each th step. The thinning factor reduces the autocorrelation of the random points in the Gibbs sampling process. As natural, this value must be an integer greater than or equal to 1.

## 4 Moments of Truncated Multivariate Elliptical Distributions

This section describes an algorithm to compute the first two moments and the variance-covariance matrix of a random vector, whose distribution belongs to the elliptical family. Furthermore, we are going to apply this algorithm to some well-known distributions. Let X be a -variate random vector that follows a truncated multivariate elliptical distribution with location parameter , positive-definite scale matrix , dgf , and support , i.e., . The more straightforward approach for this problem is to use Monte Carlo integration. Following this approach, the estimates are given by

 ˆE(X)=1nn∑i=1xi,ˆE(XX⊤)=1nn∑i=1xix⊤i,ˆCov(X)=ˆE(XX⊤)−ˆE(% X)ˆE(X)⊤, (5)

where is the th sample of the random vector X draws from . However, it is well-known that the execution time needed to perform Monte Carlo integration depends on the algorithm employed to draw samples, the number of random points () used in the approximation, and the length of the random vector (). Then, it depends on some variables that might represent a considerable computational effort. Nevertheless, we can save time when the random vector X has non-truncated components following the idea of galarza2020moments. They proposed to decompose X into two vectors, and , in such a way that is the random vector of truncated variables and is the non-truncated part, and then compute the moments for the truncated variables using any method and the remaining moments using properties of the conditional expectation. Before showing our algorithm, we state an extremely important result.

###### Proposition 4.1 (Marginal and conditional distribution of the Elliptical family)

Let be partitioned into two vectors, and , such that and has joint multivariate elliptical distribution as follows

 X=(X1X2)∼Eℓp1+p2(μ=(μ1μ2),Σ=([]ccΣ11Σ12Σ21Σ22);g(p1+p2)),

where , are location vectors, are dispersion matrices, and is the dgf. fang2018symmetric demonstrated that the elliptical family of distributions is closed under marginalization and conditioning. Hence, the distribution of and are also elliptical, with

Therefore, considering that is the vector of truncated variables with truncation region and is the vector of non-truncated variables, by Proposition 4.1 we have that

 X1∼TEℓp1(μ1,Σ11;g(p1)1,A1)and X2|(X1=x)∼Eℓp2(μ2+Σ21Σ−111(x−μ1),Σ22−Σ21Σ−111Σ12;g(p2)x).

Let and . Then, it follows that , that is

 (6)

On the other hand, we have that , with

• [leftmargin=0.4cm]

where is the expected value of a function of depending on the conditional dgf . So, the variance-covariance matrix of X is given by

 Cov(X|X∈A)=⎛⎝Ω11Ω11Σ−111Σ12Σ21Σ−111Ω11ω2.1Σ22−Σ21Σ−111(ω2.1Ip1−Ω11Σ−111)Σ12⎞⎠, (7)

Thereby, we just need Monte Carlo integration to approximate , , and (if necessary). A brief summary of how our algorithm works is given in Algorithm 2.

### 4.1 Mean and Variance for the Truncated Elliptical Distributions

Now, in this subsection, we analyze how Algorithm 2 works for some specific distributions considering all the conditions used previously.

• [leftmargin=0.4cm]

• Normal: If , the marginal distribution is and the conditional distribution is , with and . Then, With the above conditions, Algorithm 2 firstly sample from the truncated multivariate normal distribution with mean , covariance matrix , truncation region and equal to 1.

• Student-: If , the marginal and conditional distributions are and , respectively, such that , , and . For this distribution , if and , if . Therefore, the algorithm samples from the truncated distribution with location parameter , scale matrix , degrees of freedom, truncation region , and computed by

 ω2.1=ν+E(δ1(X1)|X1∈A1)ν+p1−2,

with . It is worth mention that for doubly truncated variables, the mean and the variance exist for all . Then, if X has at least two doubly truncated variables, the mean and the variance-covariance matrix exist for all . For more details about the existences of the moments see galarza2020moments.

• Pearson VII: If , then and . In this case, , if and , if . The marginal and the conditional distributions are and , respectively, such that , and . So, the proposed algorithm was implemented by sampling from the truncated multivariate Pearson VII distribution with location parameter , scale matrix , additional parameters , , and truncation region . The constant is

 ω2.1=ν+E(δ1(X1)|X1∈A1)2m−p2−2,

where is given as in the Student- distribution. For this distribution, first and second moments for doubly truncated variables exist for all . Then, if X has at least two doubly truncated variables, the mean and the variance exist for all . For more details about the existence of the moments, see Appendix B.

• Slash: If , then and . In this case, , if . The marginal distribution is and the conditional distribution is , such that , , and . Note that does not follow slash distribution, but its distribution belongs to the elliptical family (see Appendix C). So, is sampled from the truncated multivariate slash distribution with location parameter , scale matrix , degrees of freedom and truncation region . The constant is given by

 ω2.1=νν−1E(SLp1(% X1;μ1,Σ11,ν−1)SLp1(X1;μ1,Σ11,ν)∣∣X1∈A1).

This constant can be also approximated via Monte Carlo integration.

• Contaminated Normal: If , then the distributions of and are and , respectively, such that , ,