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 ExpectationMaximization (EM)
(dempster1977maximum) are employed frequently in multivariate censored data analysis under a likelihoodbased 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 mixedeffects 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 higherorder moments. On the other hand, for the TMVT distribution, the packages TTmoment (ho2015r) and MomTrunccompute its two first moments. Moreover, the first library only handles integer degrees of freedom greater than 4, while the latter can compute even highorder 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 nontruncated 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).
neal2003sliceproposed 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 nontruncated 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 skewnormal, 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 , positivedefinite scale matrix , and density generating function , if its pdf is given by
(1) 
where is a nonnegative Lebesgue measurable function on such that and denotes the determinant of matrix . Moreover,
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 variancecovariance matrix , arises when the dgf takes the form .

The multivariate Studentt 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 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
(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 positivedefinite (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 nonnegative 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 MetropolisHastings 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, positivedefinite 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
(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
(4) 
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:
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 .

[leftmargin=0.6cm]

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

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

Combining items 1 and 2, we obtain that , where

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.
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 positivedefinite 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 nonnegative 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 NewtonRaphson 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 closedform 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.
Another distribution that belongs to the elliptical family is the Kotztype 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 Kotztype distribution with settings as before, and extra parameters , and .
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 variancecovariance matrix of a random vector, whose distribution belongs to the elliptical family. Furthermore, we are going to apply this algorithm to some wellknown distributions. Let X be a variate random vector that follows a truncated multivariate elliptical distribution with location parameter , positivedefinite 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
(5) 
where is the th sample of the random vector X draws from . However, it is wellknown 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 nontruncated 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 nontruncated 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
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 nontruncated variables, by Proposition 4.1 we have that
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 variancecovariance matrix of X is given by
(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
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 variancecovariance 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
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
This constant can be also approximated via Monte Carlo integration.

Contaminated Normal: If , then the distributions of and are and , respectively, such that , ,
Comments
There are no comments yet.