# Moments and Random Number Generation for the Truncated Elliptical Family of
Distributions

Katherine A. L. Valeriano
Departamento de Estatística
Universidade Estadual de Campinas
Campinas, Brazil 13083-970
<EMAIL_ADDRESS>
&Christian E. Galarza
Departamento de Matemáticas
Escuela Superior Politécnica del Litoral
Guayaquil, Ecuador 090112
<EMAIL_ADDRESS>
&Larissa A. Matos
Departamento de Estatística
Universidade Estadual de Campinas
Campinas, Brazil 13083-970
<EMAIL_ADDRESS>

###### Abstract

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

_K_ eywords Elliptical distributions $\cdot$ Slice sampling algorithm $\cdot$
Truncated distributions $\cdot$ Truncated moments

##  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) (Dempster et al., 1977) 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, Matos et al. (2013) and Matos et al.
(2016) estimated 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 $t$ (TMVT) and a truncated
multivariate normal (TMVN) distributions, respectively.

In this context, there are a few libraries in R (R Core Team, 2021) which
provide truncated multivariate moments. For instance, the package tmvtnorm
(Wilhelm, 2015) computes the mean and the variance of the TMVN distribution by
deriving its moment generating function, while the MomTrunc library (Galarza
et al., 2021) uses a recursive approach method proposed by Kan and Robotti
(2017) to compute arbitrary higher-order moments. On the other hand, for the
TMVT distribution, the packages TTmoment (Ho et al., 2015) 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 (Galarza et al., 2021).

Variations of the EM algorithm such as Stochastic Approximation EM (SAEM)
(Delyon et al., 1999) and Monte Carlo EM (MCEM) (Wei and Tanner, 1990) replace
the conditional expectations by an approximation that requires to draw
independent random observations from a truncated distribution. For instance,
Lachos et al. (2017) 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, Lachos et al. (2019)
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, Gelfand et al. (1992) 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).
Neal (2003) 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-$t$,
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-$t$ 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 Galarza et al., 2021; de Alencar et al., 2021).

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 Muirhead (2009) and Fang (2018), a random vector
$\textbf{X}\in\mathbb{R}^{p}$ is said to follow an elliptical distribution
with location parameter $\boldsymbol{\mu}\in\mathbb{R}^{p}$, positive-definite
scale matrix $\boldsymbol{\Sigma}\in\mathbb{R}^{p\times p}$, and density
generating function $g$, if its pdf is given by

$f_{\textbf{X}}(\textbf{x})=c_{p}|\boldsymbol{\Sigma}|^{-1/2}g\left((\textbf{x}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\textbf{x}-\boldsymbol{\mu})\right),\quad\textbf{x}\in\mathbb{R}^{p},$
(1)

where $g(t)$ is a non-negative Lebesgue measurable function on $[0,\infty)$
such that $\int_{0}^{\infty}t^{p/2-1}g(t)dt<\infty$ and
$|\boldsymbol{\Sigma}|$ denotes the determinant of matrix
$\boldsymbol{\Sigma}$. Moreover,

$c_{p}=\frac{\Gamma(p/2)}{\pi^{p/2}}\left(\int_{0}^{\infty}t^{p/2-1}g(t)dt\right)^{-1}$

is the normalizing constant, with $\Gamma(\cdot)$ representing the complete
gamma function. We will use the notation $\textbf{X}\sim
E\ell_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma};g)$.

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

  * •

The multivariate normal distribution,
$\textbf{X}\sim\mathcal{N}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma})$, with
mean $\boldsymbol{\mu}$ and variance-covariance matrix $\boldsymbol{\Sigma}$,
arises when the dgf takes the form $g(t)=\exp(-t/2),t\geq 0$.

  * •

The multivariate Student-t distribution, $\textbf{X}\sim
t_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},\nu)$, where $\boldsymbol{\mu}$ is
the location parameter, $\boldsymbol{\Sigma}$ is the scale matrix, and $\nu$
is called the degrees of freedom, is obtained when
$g(t)=(1+t/\nu)^{-(\nu+p)/2},t\geq 0$.

  * •

The multivariate power exponential,
$\textbf{X}\sim\mbox{PE}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},\beta)$,
with kurtosis parameter $\beta>0$. In this case,
$g(t)=\exp(-t^{\beta}/2),t\geq 0$. A particular case of the power exponential
distribution is the normal distribution, which arises when $\beta=1$.

  * •

The multivariate slash,
$\textbf{X}\sim\mbox{SL}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},\nu)$, we
get a random variable with multivariate slash distribution when
$g(t)=\int_{0}^{1}u^{\nu+p/2-1}\exp\\{-ut/2\\}du,t\geq 0,\nu>0$.

  * •

The multivariate Pearson VII distribution,
$\textbf{X}\sim\mbox{PVII}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},m,\nu)$,
with parameters $\boldsymbol{\mu}\in\mathbb{R}^{p}$,
$\boldsymbol{\Sigma}\in\mathbb{R}^{p\times p}$, $m>p/2$, and $\nu>0$ is
obtained when $g(t)=(1+t/\nu)^{-m},t\geq 0$.

For more distributions belonging to this family, please see Fang (2018).

###  2.2 Truncated Elliptical Family of Distributions

Let ${A}\subseteq\mathbb{R}^{p}$ be a measurable set. We say that a random
vector $\textbf{Y}\in\mathbb{R}^{p}$ has truncated elliptical distribution
with support ${A}$, location parameter $\boldsymbol{\mu}\in\mathbb{R}^{p}$,
scale parameter $\boldsymbol{\Sigma}\in\mathbb{R}^{p\times p}$ and dgf $g$, if
its pdf is given by

$f_{\textbf{Y}}(\textbf{y})=\frac{g\left((\textbf{y}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\textbf{y}-\boldsymbol{\mu})\right)}{\int_{A}g\left((\textbf{y}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\textbf{y}-\boldsymbol{\mu})\right)d\textbf{y}}=\frac{f_{\textbf{X}}(\textbf{y})}{\mbox{Pr}(\textbf{X}\in{A})},\quad\textbf{y}\in{A},$
(2)

where
$\textbf{X}\sim\mbox{E}\ell_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma};g)$. We
use the notation
$\textbf{Y}\sim\mbox{TE}\ell(\boldsymbol{\mu},\boldsymbol{\Sigma};g,A)$.
Notice that the pdf of Y is written as the ratio between the pdf of
$\textbf{X}\sim\mbox{E}\ell_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma};g)$ and
$\mbox{Pr}(\textbf{X}\in{A})$, so the pdf of Y exists if the pdf of X does,
which occurs if $\boldsymbol{\Sigma}$ is a positive-definite (see, Morán-
Vásquez and Ferrari, 2019). The variable Y is also said to be an elliptical
distribution truncated on $A$, being represented by
$\textbf{Y}=\textbf{X}|\,(\textbf{X}\in A)$.

As in the elliptical family of distributions, the dgf $g$ determines any
distribution within the truncated elliptical class of distributions, for
example, if $g(t)=(1+t/\nu)^{-(\nu+p)/2},t\geq 0,\nu>0$, then Y has TMVT
distribution. We will denote the different members of the truncated elliptical
family defined in the subsection before as
$\textbf{Y}\sim\mbox{TN}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma};A)$ for the
TMVN distribution,
$\textbf{Y}\sim\mbox{T}t_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},\nu;A)$ for
the TMVT distribution,
$\textbf{Y}\sim\mbox{TPE}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},\beta;A)$
for the truncated multivariate power exponential,
$\textbf{Y}\sim\mbox{TSL}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},\nu;A)$ for
the truncated multivariate slash distribution, and
$\textbf{Y}\sim\mbox{TPVII}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},m,\nu;A)$
for the truncated multivariate Pearson VII distribution.

###  2.3 Slice Sampling Algorithm

Introduced by Neal (2003), 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
$(p+1)$-dimensional region under the graph of $f(\textbf{x})$, a non-negative
function proportional to the pdf of X. Hence, let $Y$ be an auxiliary variable
such that the joint pdf of X and $Y$ is uniform over the region
$U=\\{(\textbf{x},y):0<y<f(\textbf{x})\\}$, i.e.,
$f_{X,Y}(\textbf{x},y)\propto\mathbb{I}\left(0<y<f(\textbf{x})\right)$, with
$\mathbb{I}(\cdot)$ being the indicator function. Therefore, we can obtain
samples from the distribution of X by sampling jointly $(\textbf{x},y)$ and
then ignoring $y$ values.

Figure 1: Slice sampling algorithm for univariate random variables.

Note that generating independent random points uniformly distributed on $U$
may not be easy. To overcome this problem, Neal (2003) 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 $k$ works as follows: given
the current value of $\textbf{x}_{k-1}$ sample $y_{k}$ from
$Y|(\textbf{X}=\textbf{x}_{k-1})\sim\mathcal{U}\left(0,f(\textbf{x}_{k-1})\right)$,
then draw $\textbf{x}_{k}$ from the conditional distribution of X given
$y_{k}$, which is uniform over the region
$S_{k}=\\{\textbf{x}:y_{k}<f(\textbf{x})\\}$, i.e.,
$\textbf{X}|(Y=y_{k})\sim\mathcal{U}(\\{\textbf{x}:y_{k}<f(\textbf{x})\\})$,
for all $k=1,2,\ldots,n$, where $n$ is the desired sample size.

Figure 1 shows the steps of the slice sampling algorithm for $X$ being a
univariate random variable. Given an initial value $X=x_{0}$, we draw $y_{1}$
uniformly over the interval $\left(0,f(x_{0})\right)$ and then we sample
$x_{1}$ from the conditional distribution of $X|(Y=y_{1})$, i.e., uniformly
over the interval $S_{1}=\\{x:y_{1}<f(x)\\}$. These two steps are repeated $n$
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
$p$-variate truncated elliptical distribution with zero location parameter,
positive-definite scale matrix $\mathbf{R}\in\mathbb{R}^{p\times p}$, dgf $g$,
and truncation region ${A}=\\{\textbf{x}:\mathbf{a}<\textbf{x}<\mathbf{b}\\}$,
$\mathbf{a}<\mathbf{b}\in\mathbb{R}^{p}$, in other words, we will consider
$\textbf{X}\sim\mbox{TE}\ell_{p}(\mathbf{0},\mathbf{R};g,A)$. Here
$\mathbf{R}$ is a correlation matrix, such that the scale matrix can be
written as
$\boldsymbol{\Sigma}=\boldsymbol{\Lambda}\mathbf{R}\boldsymbol{\Lambda}$,
where $\boldsymbol{\Lambda}=\mathrm{diag}(\sigma_{11},\ldots,\sigma_{pp})$.
The pdf of X is given by

$f_{\textbf{X}}(\textbf{x})\propto
g\left(\textbf{x}^{\top}\textbf{R}^{-1}\textbf{x}\right)\mathbb{I}\left(\textbf{x}\in{A}\right),$
(3)

Now, in order to sample uniformly from the $(p+1)$-dimensional region under
the plot of $f_{\textbf{X}}(\textbf{x})$, we introduce an auxiliary variable
$Y$, such that the joint pdf of X and $Y$ is

$f_{\textbf{X},Y}(\textbf{x},y)\propto\mathbb{I}\left(0<y<g\left(\textbf{x}^{\top}\textbf{R}^{-1}\textbf{x}\right)\right)\mathbb{I}\left(\textbf{a}<\textbf{x}<\textbf{b}\right).$
(4)

It is enough to calculate the conditional distributions of $Y|\textbf{X}$ and
$\textbf{X}|Y$ 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:

$\displaystyle f_{Y|\textbf{X}}(y|\textbf{x})$ $\displaystyle\propto$
$\displaystyle\mathbb{I}\left(0<y<g\left(\textbf{x}^{\top}\textbf{R}^{-1}\textbf{x}\right)\right)\quad\mbox{
and }\quad$ $\displaystyle f_{\textbf{X}|Y}(\textbf{x}|y)$
$\displaystyle\propto$
$\displaystyle\mathbb{I}\left(\big{\\{}\textbf{x}:y<g\left(\textbf{x}^{\top}\textbf{R}^{-1}\textbf{x}\right)\big{\\}}\cap\\{\textbf{a}<\textbf{x}<\textbf{b}\\}\right).$

Note that sampling $y$ from the distribution of $Y|(\textbf{X}=\textbf{x})$ is
straightforward, but sampling from $\textbf{X}|(Y=y)$ is not trivial. Thus, we
use the idea of Ho et al. (2012), that consists in sampling each element of X
given the remaining elements, i.e., sampling $X_{j}$ given
$\textbf{x}_{-j}=(x_{1},\ldots,x_{j-1},x_{j+1},\ldots,x_{p})^{\top}$ and $y$,
for all $j=1,\ldots,p$. Hence, the following steps are performed to draw a
random number from the distribution of $X_{j}|\textbf{x}_{-j},y$.

  1. 1.

Let $\kappa_{y}=g^{-1}(y)$. Since $g$ is a strictly decreasing function, it
follows that $y<g(\textbf{x}^{\top}\textbf{R}^{-1}\textbf{x})$ is equivalent
to $\kappa_{y}>\textbf{x}^{\top}\textbf{R}^{-1}\textbf{x}$.

  2. 2.

Write
$\textbf{x}^{\top}\textbf{R}^{-1}\textbf{x}=\rho^{jj}\left(x_{j}-\lambda_{j}\right)^{2}-\rho^{jj}\lambda_{j}^{2}+\eta_{j}$,
where $\rho^{ij}$ is the $(i,j)$th element of the inverse of R,
$\eta_{j}=\sum_{t\neq j}\sum_{r\neq j}x_{t}x_{r}\rho^{tr}$ and
$\lambda_{j}=-\frac{1}{\rho^{jj}}\sum_{r\neq j}x_{r}\rho^{jr}$.

  3. 3.

Combining items 1 and 2, we obtain that
$\lambda_{j}-\tau_{j}<x_{j}<\lambda_{j}+\tau_{j}$, where
$\tau_{j}=\left(\lambda_{j}^{2}+\frac{1}{\rho^{jj}}\left(\kappa_{y}-\eta_{j}\right)\right)^{1/2}$

  4. 4.

Because $x_{j}\in(a_{j},b_{j})$, thereby
$a_{j}^{*}=\max(a_{j},\lambda_{j}-\tau_{j})<x_{j}<\min(b_{j},\lambda_{j}+\tau_{j})=b_{j}^{*}$.

Therefore, the steps to draw $n$ samples from a $p$-variate truncated
elliptical distribution
$\textbf{X}\sim\mbox{TE}\ell_{p}(\textbf{0},\textbf{R};g,A)$ 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 $g$ is strictly decreasing has been used in step 1. A general case can be
easily considered by studying the extrema points of $g$.

Input: Sample size $n\geq 1$, initial value $\textbf{x}_{0}\in\mathbb{R}^{p}$,
scale matrix $\textbf{R}\in\mathbb{R}^{p\times p}$, lower bound
$\textbf{a}\in\mathbb{R}^{p}$, upper bound $\textbf{b}\in\mathbb{R}^{p}$ and
strictly decreasing dgf $g(t),t\geq 0$.

Initialization;

for  _$i\leftarrow 1$ to $n$_ do

Sample $y_{i}$ from
$Y|\textbf{x}_{i-1}\sim\mathcal{U}(0,g(\textbf{x}_{i-1}^{\top}\textbf{R}^{-1}\textbf{x}_{i-1}))$;

$\kappa_{y}\leftarrow g^{-1}(y_{i})$;

for  _$j\leftarrow 1$ to $p$_ do

$\eta_{j}\leftarrow\displaystyle\sum_{t\neq j}\sum_{r\neq
j}x_{t}x_{r}\rho^{tr};\quad$
$\lambda_{j}\leftarrow\displaystyle-\frac{1}{\rho^{jj}}\sum_{r\neq
j}x_{r}\rho^{jr};\quad$
$\tau_{j}\leftarrow\displaystyle\left(\lambda_{j}^{2}+\frac{1}{\rho^{jj}}\left(\kappa_{y}-\eta_{j}\right)\right)^{1/2}$;

Draw $x_{j}$ from
$X_{j}|\textbf{x}_{-j},y\sim\mathcal{U}(\max(a_{j},\lambda_{j}-\tau_{j}),\min(b_{j},\lambda_{j}+\tau_{j}))$;

$\textbf{x}_{i}[j]\leftarrow x_{j};\quad$ $\textbf{X}[i,j]\leftarrow x_{j}$;

end for

end for

Result: X

Algorithm 1 Slice sampling algorithm

Moreover, members of the truncated elliptical family of distributions are
closed under affine transformations (Fang, 2018). Hence drawing samples from
$\textbf{Y}\sim\mbox{TE}\ell_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma};g,(\textbf{a},\textbf{b}))$
may be readily done by sampling first from
$\textbf{X}\sim\mbox{TE}\ell_{p}(\textbf{0},\textbf{R};g,(\textbf{a}^{*},\textbf{b}^{*}))$,
and then recovery Y by the following transformation
$\textbf{Y}=\boldsymbol{\mu}+\boldsymbol{\Lambda}\textbf{X}$, such that
$\boldsymbol{\Sigma}=\boldsymbol{\Lambda}\textbf{R}\boldsymbol{\Lambda}$,
$\textbf{a}^{*}=\boldsymbol{\Lambda}^{-1}\left(\textbf{a}-\boldsymbol{\mu}\right)$,
and
$\textbf{b}^{*}=\boldsymbol{\Lambda}^{-1}\left(\textbf{b}-\boldsymbol{\mu}\right)$.

###  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, $n\geq 1$ 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-$t$, 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
$\boldsymbol{\mu}=(0,0)^{\top}$, scale matrix elements
$\sigma_{11}=\sigma_{22}=1$, and $\sigma_{12}=\sigma_{21}=0.70$, and
truncation region ${A}=\\{\textbf{x}:\textbf{a}<\textbf{x}<\textbf{b}\\}$,
with $\textbf{a}=(-2,-2)^{\top}$ and $\textbf{b}=(3,2)^{\top}$. The
distributions considered are the predefined ones in the package.

  * •

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-$t$ with $\nu=3$ 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 $\beta=2$

⬇

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 $m=5/2$ and $\nu=3$

⬇

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 $\nu=0.70$ and $\rho=0.20$

⬇

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 $\kappa_{y}=g^{-1}(y)$ for the
normal, Student-$t$, 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
$g(t)-y=0,t\geq 0$, through the Newton-Raphson algorithm for the contaminated
normal, and using Brent’s method (Brent, 2013), 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 (Andersen et al., 2020), and the R
environment. For instance, for the dgf $g(t)=e^{-t}$, 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 $g(t)$, 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
$g(t)=e^{-t}/(1+e^{-t})^{2},t\geq 0$, 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 $r>0,s>0$, and $2N+p>2$, whose dgf is
$g(t)=t^{N-1}e^{-rt^{s}},t\geq 0$ (Fang, 2018). For this distribution, $g(t)$
is not strictly decreasing, however, for $(2-p)/2<N\leq 1$, it holds. Hence,
our proposal works for $r>0$, $s>0$, and $(2-p)/2<N\leq 1$. 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 $r=2,s=1/4$, and $N=1/2$.

⬇

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 $n=10^{4}$
observations sampled from each of the truncated bivariate distributions
referred above.

Figure 2: Scatterplot and marginal histograms for the $n=10^{4}$ observations
sampled for some bivariate truncated elliptical distributions.

As mentioned by Robert and Casella (2010) and Ho et al. (2012), 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 $k$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 $p$-variate random vector that
follows a truncated multivariate elliptical distribution with location
parameter $\boldsymbol{\mu}\in\mathbb{R}^{p}$, positive-definite scale matrix
$\boldsymbol{\Sigma}\in\mathbb{R}^{p\times p}$, dgf $g$, and support
${A}\subseteq\mathbb{R}^{p}$, i.e.,
$\textbf{X}\sim\mbox{TE}\ell_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma};g,A)$.
The more straightforward approach for this problem is to use Monte Carlo
integration. Following this approach, the estimates are given by

$\widehat{\mathbb{E}(\textbf{X})}=\frac{1}{n}\sum_{i=1}^{n}\textbf{x}_{i},\quad\widehat{\mathbb{E}(\textbf{X}\textbf{X}^{\top})}=\frac{1}{n}\sum_{i=1}^{n}\textbf{x}_{i}\textbf{x}_{i}^{\top},\quad\widehat{\mathrm{Cov}(\textbf{X})}=\widehat{\mathbb{E}(\textbf{X}\textbf{X}^{\top})}-\widehat{\mathbb{E}(\textbf{X})}\widehat{\mathbb{E}(\textbf{X})}^{\top},$
(5)

where $\textbf{x}_{i}$ is the $i$th sample of the random vector X draws from
$\mbox{TE}\ell_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma};g,A)$. 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
($n$) used in the approximation, and the length of the random vector ($p$).
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 Galarza et al. (2020). They
proposed to decompose X into two vectors, $\textbf{X}_{1}$ and
$\textbf{X}_{2}$, in such a way that $\textbf{X}_{1}$ is the random vector of
truncated variables and $\textbf{X}_{2}$ 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 $\textbf{X}\in\mathbb{R}^{p}$ be partitioned into two vectors,
$\textbf{X}_{1}\in\mathbb{R}^{p_{1}}$ and
$\textbf{X}_{2}\in\mathbb{R}^{p_{2}}$, such that $p=p_{1}+p_{2}$ and
$\textbf{X}=(\textbf{X}_{1}^{\top},\textbf{X}_{2}^{\top})^{\top}$ has joint
multivariate elliptical distribution as follows

$\textbf{X}=\left(\begin{array}[]{c}\textbf{X}_{1}\\\
\textbf{X}_{2}\end{array}\right)\sim\mbox{E}\ell_{p_{1}+p_{2}}\left(\boldsymbol{\mu}=\left(\begin{array}[]{c}\boldsymbol{\mu}_{1}\\\
\boldsymbol{\mu}_{2}\end{array}\right),\boldsymbol{\Sigma}=\left(\begin{array}[]{cc}\boldsymbol{\Sigma}_{11}&\boldsymbol{\Sigma}_{12}\\\
\boldsymbol{\Sigma}_{21}&\boldsymbol{\Sigma}_{22}\end{array}\right);g^{(p_{1}+p_{2})}\right),$

where $\boldsymbol{\mu}_{1}\in\mathbb{R}^{p_{1}}$,
$\boldsymbol{\mu}_{2}\in\mathbb{R}^{p_{2}}$ are location vectors,
$\boldsymbol{\Sigma}_{11}\in\mathbb{R}^{p_{1}\times
p_{1}},\boldsymbol{\Sigma}_{22}\in\mathbb{R}^{p_{2}\times
p_{2}},\boldsymbol{\Sigma}_{12}\in\mathbb{R}^{p_{1}\times
p_{2}},\boldsymbol{\Sigma}_{21}\in\mathbb{R}^{p_{2}\times p_{1}}$ are
dispersion matrices, and $g^{(p_{1}+p_{2})}$ is the dgf. Fang (2018)
demonstrated that the elliptical family of distributions is closed under
marginalization and conditioning. Hence, the distribution of $\textbf{X}_{1}$
and $\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x})$ are also elliptical, with
$\textbf{X}_{1}\sim\mbox{E}\ell_{p_{1}}\left(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11};g_{1}^{(p_{1})}\right),$
$\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x})\sim\mbox{E}\ell_{p_{2}}\left(\boldsymbol{\mu}_{2}+\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}(\textbf{x}-\boldsymbol{\mu}_{1}),\boldsymbol{\Sigma}_{22}-\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12};g^{(p_{2})}_{\textbf{x}}\right).$

Therefore, considering that $\textbf{X}_{1}$ is the vector of truncated
variables with truncation region ${A}_{1}$ and $\textbf{X}_{2}$ is the vector
of non-truncated variables, by Proposition 4.1 we have that

$\displaystyle\textbf{X}_{1}\sim\mbox{TE}\ell_{p_{1}}(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11};g_{1}^{(p_{1})},A_{1})\quad\mbox{and}$
$\displaystyle\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x})\sim\mbox{E}\ell_{p_{2}}\left(\boldsymbol{\mu}_{2}+\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}(\textbf{x}-\boldsymbol{\mu}_{1}),\boldsymbol{\Sigma}_{22}-\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12};g^{(p_{2})}_{\textbf{x}}\right).$

Let $\boldsymbol{\xi}_{1}=\mathbb{E}(\textbf{X}_{1}|\textbf{X}_{1}\in{A}_{1})$
and
$\boldsymbol{\Omega}_{11}=\mathrm{Cov}(\textbf{X}_{1}|\textbf{X}_{1}\in{A}_{1})$.
Then, it follows that
$\mathbb{E}(\textbf{X}|\textbf{X}\in{A})=\mathbb{E}(\mathbb{E}(\textbf{X}|\textbf{X}_{1})|\textbf{X}_{1}\in{A}_{1})$,
that is

$\displaystyle\mathbb{E}(\textbf{X}|\textbf{X}\in{A})=\mathbb{E}\left(\begin{array}[]{c}\textbf{X}_{1}\\\
\boldsymbol{\mu}_{2}+\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}(\textbf{X}_{1}-\boldsymbol{\mu}_{1})\end{array}\Big{|}\textbf{X}_{1}\in{A}_{1}\right)=\left(\begin{array}[]{c}\boldsymbol{\xi}_{1}\\\
\boldsymbol{\mu}_{2}+\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}(\boldsymbol{\xi}_{1}-\boldsymbol{\mu}_{1})\end{array}\right).$
(10)

On the other hand, we have that
$\mathrm{Cov}(\textbf{X}|\textbf{X}\in{A})=\mathrm{Cov}(\mathbb{E}(\textbf{X}|\textbf{X}_{1})|\textbf{X}_{1}\in{A}_{1})+\mathbb{E}(\mathrm{Cov}(\textbf{X}|\textbf{X}_{1})|\textbf{X}_{1}\in{A}_{1})$,
with

  * •

$\mathrm{Cov}(\textbf{X}_{1},\mathbb{E}(\textbf{X}_{2}|\textbf{X}_{1})|\textbf{X}_{1}\in{A}_{1})=\mathrm{Cov}(\textbf{X}_{1},\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\textbf{X}_{1}|\textbf{X}_{1}\in{A}_{1})=\boldsymbol{\Omega}_{11}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12},$

  * •

$\mathrm{Cov}(\mathbb{E}(\textbf{X}_{2}|\textbf{X}_{1})|\textbf{X}_{1}\in{A}_{1})=\mathrm{Cov}(\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\textbf{X}_{1}|\textbf{X}_{1}\in{A}_{1})=\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Omega}_{11}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12},$

  * •

$\mathbb{E}(\mathrm{Cov}(\textbf{X}_{2}|\textbf{X}_{1})|\textbf{X}_{1}\in{A}_{1})=\omega_{2.1}(\boldsymbol{\Sigma}_{22}-\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12}),$

where
$\omega_{2.1}=\mathbb{E}\left(h(\textbf{X}_{1})|\textbf{X}_{1}\in{A}_{1}\right)$
is the expected value of a function $h$ of $\textbf{X}_{1}$ depending on the
conditional dgf $g_{\textbf{x}_{1}}^{(p_{2})}$. So, the variance-covariance
matrix of X is given by

$\displaystyle\mathrm{Cov}(\textbf{X}|\textbf{X}\in{A})=\left(\begin{array}[]{cc}\boldsymbol{\Omega}_{11}&\boldsymbol{\Omega}_{11}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12}\\\
\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Omega}_{11}&\omega_{2.1}\boldsymbol{\Sigma}_{22}-\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\left(\omega_{2.1}\textbf{I}_{p_{1}}-\boldsymbol{\Omega}_{11}\boldsymbol{\Sigma}_{11}^{-1}\right)\boldsymbol{\Sigma}_{12}\end{array}\right),$
(13)

Thereby, we just need Monte Carlo integration to approximate
$\boldsymbol{\xi}_{1}$, $\boldsymbol{\Omega}_{11}$, and $\omega_{2.1}$ (if
necessary). A brief summary of how our algorithm works is given in Algorithm
2.

Input: Sample size $n\geq 1$, location parameter
$\boldsymbol{\mu}\in\mathbb{R}^{p}$, scale matrix
$\boldsymbol{\Sigma}\in\mathbb{R}^{p\times p}$, lower bound
$\textbf{a}\in\mathbb{R}^{p}$, upper bound $\textbf{b}\in\mathbb{R}^{p}$ and
dgf $g(t),t\geq 0$.

Identify:
$\boldsymbol{\mu}_{1},\boldsymbol{\mu}_{2},\boldsymbol{\Sigma}_{11},\boldsymbol{\Sigma}_{22},\boldsymbol{\Sigma}_{12},{A}_{1}=\\{\textbf{x}_{1}:\textbf{a}_{1}<\textbf{x}_{1}<\textbf{b}_{1}\\}$;

Draw $\textbf{x}_{1i}$ from
$\textbf{X}_{1}\sim\mbox{TE}\ell_{p_{1}}(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11};g^{p_{1}},A_{1}),\quad
i=1,\ldots,n$;

$\widehat{\boldsymbol{\xi}}_{1}\leftarrow\displaystyle\frac{1}{n}\sum_{i=1}^{n}\textbf{x}_{1i};\quad$
$\widehat{\boldsymbol{\Omega}}_{11}\leftarrow\displaystyle\frac{1}{n}\sum_{i=1}^{n}\textbf{x}_{1i}\textbf{x}_{1i}^{\top}-\widehat{\boldsymbol{\xi}}_{1}\widehat{\boldsymbol{\xi}}_{1}^{\top};\quad$
$\widehat{\omega}_{2.1}\leftarrow\displaystyle\frac{1}{n}\sum_{i=1}^{n}h(\textbf{x}_{1i})$;

$\widehat{\mathbb{E}(\textbf{X})}\leftarrow\left(\begin{array}[]{cc}\widehat{\boldsymbol{\xi}}_{1}\\\
\boldsymbol{\mu}_{2}+\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}(\widehat{\boldsymbol{\xi}}_{1}-\boldsymbol{\mu}_{1})\end{array}\right)$;

$\widehat{\mathrm{Cov}(\textbf{X})}\leftarrow\left(\begin{array}[]{cc}\widehat{\boldsymbol{\Omega}}_{11}&\widehat{\boldsymbol{\Omega}}_{11}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12}\\\
\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\widehat{\boldsymbol{\Omega}}_{11}&\widehat{\omega}_{2.1}\boldsymbol{\Sigma}_{22}-\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\left(\widehat{\omega}_{2.1}\textbf{I}_{p_{1}}-\widehat{\boldsymbol{\Omega}}_{11}\boldsymbol{\Sigma}_{11}^{-1}\right)\boldsymbol{\Sigma}_{12}\end{array}\right)$;

$\widehat{\mathbb{E}(\textbf{X}\textbf{X}^{\top})}\leftarrow\widehat{\mathrm{Cov}(\textbf{X})}+\widehat{\mathbb{E}(\textbf{X})}\widehat{\mathbb{E}(\textbf{X})}^{\top}$;

Result:
$\widehat{\mathbb{E}(\textbf{X})},\widehat{\mathbb{E}(\textbf{X}\textbf{X}^{\top})},\widehat{\mathrm{Cov}(\textbf{X})}$

Algorithm 2 Mean and variance approximation

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

  * •

Normal: If
$\textbf{X}\sim\mathcal{N}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma})$, the
marginal distribution is
$\textbf{X}_{1}\sim\mathcal{N}_{p_{1}}(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11})$
and the conditional distribution is
$\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x})\sim\mathcal{N}_{p_{2}}(\boldsymbol{\mu}_{2.1},\boldsymbol{\Sigma}_{2.1})$,
with
$\boldsymbol{\mu}_{2.1}=\boldsymbol{\mu}_{2}+\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}(\textbf{x}-\boldsymbol{\mu}_{1})$
and
$\boldsymbol{\Sigma}_{2.1}=\boldsymbol{\Sigma}_{22}-\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12}$.
Then, With the above conditions, Algorithm 2 firstly sample $\textbf{X}_{1}$
from the truncated multivariate normal distribution with mean
$\boldsymbol{\mu}_{1}$, covariance matrix $\boldsymbol{\Sigma}_{11}$,
truncation region ${A}_{1}$ and $\omega_{2.1}$ equal to 1.

  * •

Student-$t$: If $\textbf{X}\sim
t_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},\nu),\nu>0$, the marginal and
conditional distributions are $\textbf{X}_{1}\sim
t_{p_{1}}(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11},\nu)$ and
$\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x})\sim
t_{p_{2}}(\boldsymbol{\mu}_{2.1},\lambda\boldsymbol{\Sigma}_{2.1},\nu+p_{1})$,
respectively, such that
$\boldsymbol{\mu}_{2.1}=\boldsymbol{\mu}_{2}+\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}(\textbf{x}-\boldsymbol{\mu}_{1})$,
$\lambda=(\nu+\delta_{1}(\textbf{x}))/(\nu+p_{1})$,
$\boldsymbol{\Sigma}_{2.1}=\boldsymbol{\Sigma}_{22}-\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12}$
and
$\delta_{1}(\textbf{x})=(\textbf{x}-\boldsymbol{\mu}_{1})^{\top}\boldsymbol{\Sigma}_{11}^{-1}(\textbf{x}-\boldsymbol{\mu}_{1})$.
For this distribution $\mathbb{E}(\textbf{X})<\infty$, if $\nu>1$ and
$\mathrm{Cov}(\textbf{X})<\infty$, if $\nu>2$. Therefore, the algorithm
samples $\textbf{X}_{1}$ from the truncated $t$ distribution with location
parameter $\boldsymbol{\mu}_{1}$, scale matrix $\boldsymbol{\Sigma}_{11}$,
$\nu$ degrees of freedom, truncation region ${A}_{1}$, and $\omega_{2.1}$
computed by

$\displaystyle\omega_{2.1}=\frac{\nu+\mathbb{E}(\delta_{1}(\textbf{X}_{1})|\textbf{X}_{1}\in{A}_{1})}{\nu+p_{1}-2},$

with
$\mathbb{E}(\delta_{1}(\textbf{X}_{1})|\textbf{X}_{1}\in{A}_{1})=\mathrm{tr}(\boldsymbol{\Omega}_{11}\boldsymbol{\Sigma}_{11}^{-1})+(\boldsymbol{\xi}_{1}-\boldsymbol{\mu}_{1})^{\top}\boldsymbol{\Sigma}_{11}^{-1}(\boldsymbol{\xi}_{1}-\boldsymbol{\mu}_{1})$.
It is worth mention that for doubly truncated variables, the mean and the
variance exist for all $\nu>0$. Then, if X has at least two doubly truncated
variables, the mean and the variance-covariance matrix exist for all $\nu>0$.
For more details about the existences of the moments see Galarza et al.
(2020).

  * •

Pearson VII: If
$\textbf{X}\sim\mbox{PVII}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},m,\nu),m>p/2,\nu>0$,
then $\mathbb{E}(\textbf{X})=\boldsymbol{\mu}$ and
$\mathrm{Cov}(\textbf{X})=\frac{\nu}{2m-p-2}\boldsymbol{\Sigma}$. In this
case, $\mathbb{E}(\textbf{X})<\infty$, if $m>(p+1)/2$ and
$\mathrm{Cov}(\textbf{X})<\infty$, if $m>(p+2)/2$. The marginal and the
conditional distributions are
$\textbf{X}_{1}\sim\mbox{PVII}_{p_{1}}(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11},m-p_{2}/2,\nu)$
and
$\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x})\sim\mbox{PVII}_{p_{2}}(\boldsymbol{\mu}_{2.1},\boldsymbol{\Sigma}_{2.1},m,\nu+\delta_{1}(\textbf{x}))$,
respectively, such that
$\boldsymbol{\mu}_{2.1}=\boldsymbol{\mu}_{2}+\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}(\textbf{x}-\boldsymbol{\mu}_{1})$,
$\boldsymbol{\Sigma}_{2.1}=\boldsymbol{\Sigma}_{22}-\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12}$
and
$\delta_{1}(\textbf{x})=(\textbf{x}-\boldsymbol{\mu}_{1})^{\top}\boldsymbol{\Sigma}_{11}^{-1}(\textbf{x}-\boldsymbol{\mu}_{1})$.
So, the proposed algorithm was implemented by sampling $\textbf{X}_{1}$ from
the truncated multivariate Pearson VII distribution with location parameter
$\boldsymbol{\mu}_{1}$, scale matrix $\boldsymbol{\Sigma}_{11}$, additional
parameters $m-p_{2}/2>p_{1}/2$, $\nu>0$, and truncation region ${A}_{1}$. The
constant $\omega_{2.1}$ is

$\omega_{2.1}=\frac{\nu+\mathbb{E}(\delta_{1}(\textbf{X}_{1})|\textbf{X}_{1}\in{A}_{1})}{2m-p_{2}-2},$

where $\mathbb{E}(\delta_{1}(\textbf{X}_{1})|\textbf{X}_{1}\in{A}_{1})$ is
given as in the Student-$t$ distribution. For this distribution, first and
second moments for doubly truncated variables exist for all $m>p/2$. Then, if
X has at least two doubly truncated variables, the mean and the variance exist
for all $m>p/2$. For more details about the existence of the moments, see
Appendix B.

  * •

Slash: If
$\textbf{X}\sim\mbox{SL}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},\nu),\nu>0$,
then $\mathbb{E}(\textbf{X})=\boldsymbol{\mu}$ and
$\mathrm{Cov}(\textbf{X})=\frac{\nu}{\nu-1}\boldsymbol{\Sigma}$. In this case,
$\mathrm{Cov}(\textbf{X})<\infty$, if $\nu>1$. The marginal distribution is
$\textbf{X}_{1}\sim\mbox{SL}_{p_{1}}(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11},\nu)$
and the conditional distribution is
$\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x})\sim\mbox{E}\ell_{p_{2}}(\boldsymbol{\mu}_{2.1},\boldsymbol{\Sigma}_{2.1};g^{(p_{2})})$,
such that
$\boldsymbol{\mu}_{2.1}=\boldsymbol{\mu}_{2}+\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}(\textbf{x}-\boldsymbol{\mu}_{1})$,
$\boldsymbol{\Sigma}_{2.1}=\boldsymbol{\Sigma}_{22}-\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12}$,
$g^{(p_{2})}(t)=\int_{0}^{1}u^{\nu+p/2-1}\exp\\{-u(t+\delta_{1}(\textbf{x}))/2\\}du$
and
$\delta_{1}(\textbf{x})=(\textbf{x}-\boldsymbol{\mu}_{1})^{\top}\boldsymbol{\Sigma}_{11}^{-1}(\textbf{x}-\boldsymbol{\mu}_{1})$.
Note that $\textbf{X}_{2}|\textbf{X}_{1}$ does not follow slash distribution,
but its distribution belongs to the elliptical family (see Appendix C). So,
$\textbf{X}_{1}$ is sampled from the truncated multivariate slash distribution
with location parameter $\boldsymbol{\mu}_{1}$, scale matrix
$\boldsymbol{\Sigma}_{11}$, $\nu$ degrees of freedom and truncation region
${A}_{1}$. The constant $\omega_{2.1}$ is given by

$\omega_{2.1}=\frac{\nu}{\nu-1}\mathbb{E}\left(\frac{\mbox{SL}_{p_{1}}(\textbf{X}_{1};\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11},\nu-1)}{\mbox{SL}_{p_{1}}(\textbf{X}_{1};\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11},\nu)}\Big{|}\textbf{X}_{1}\in{A}_{1}\right).$

This constant can be also approximated via Monte Carlo integration.

  * •

Contaminated Normal: If
$\textbf{X}\sim\mbox{CN}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},\nu,\rho),0<\nu,\rho<1$,
then the distributions of $\textbf{X}_{1}$ and
$\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x})$ are
$\mbox{CN}_{p_{1}}(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11},\nu,\rho)$
and
$\mbox{CN}_{p_{2}}(\boldsymbol{\mu}_{2.1},\boldsymbol{\Sigma}_{2.1},\nu_{2.1},\rho)$,
respectively, such that
$\boldsymbol{\mu}_{2.1}=\boldsymbol{\mu}_{2}+\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}(\textbf{x}-\boldsymbol{\mu}_{1})$,
$\boldsymbol{\Sigma}_{2.1}=\boldsymbol{\Sigma}_{22}-\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12}$,
$\nu_{2.1}=\nu\phi_{p_{1}}(\textbf{x};\boldsymbol{\mu}_{1},\rho^{-1}\boldsymbol{\Sigma}_{11})/(\nu\phi_{p_{1}}(\textbf{x};\boldsymbol{\mu}_{1},\rho^{-1}\boldsymbol{\Sigma}_{11})+(1-\nu)\phi_{p_{1}}(\textbf{x};\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11}))$
and $\phi_{p}(\textbf{x};\boldsymbol{\mu},\boldsymbol{\Sigma})$ denotes the
pdf of a $p$-variate normal distribution with mean $\boldsymbol{\mu}$,
variance matrix $\boldsymbol{\Sigma}$ evaluated at point
$\textbf{x}\in\mathbb{R}^{p_{1}}$. Thus, $\textbf{X}_{1}$ is sampled from the
truncated contaminated normal distribution with parameters
$\boldsymbol{\mu}_{1}$, $\boldsymbol{\Sigma}_{11}$, $\nu$ and $\rho$. The
constant is $\omega_{2.1}$ is
$\omega_{2.1}={\nu}_{2.1}^{*}/\rho+1-{\nu}_{2.1}^{*}$, where
${\nu}_{2.1}^{*}=\mathbb{E}(\nu_{2.1}|\textbf{X}_{1}\in{A}_{1})$, this value
is also approximated via Monte Carlo integration.

  * •

Power exponential: If
$\textbf{X}\sim\mbox{PE}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},\beta),\beta>0$,
then $\mathbb{E}(\textbf{X})=\boldsymbol{\mu}$ and
$Cov(\textbf{X})=\omega\boldsymbol{\Sigma}$, with
$\omega=2^{1/\beta}\Gamma(\frac{n+2}{2\beta})/(n\Gamma(\frac{n}{2\beta}))$.
The marginal distribution of $\textbf{X}_{1}$ belongs to the elliptical family
of distributions with dgf
$g^{(p_{1})}(t)=t^{\frac{p-p_{1}}{2}}\int_{0}^{1}w^{\frac{p_{1}-p}{2}}(1-w)^{\frac{p-p_{1}}{2}-1}\exp\\{-\frac{t^{\beta}}{2w^{\beta}}\\}dw$,
$\textbf{X}_{1}\sim\mbox{E}\ell_{p_{1}}(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11};g^{(p_{1})})$,
and the conditional distribution is
$\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x})\sim\mbox{E}\ell_{p_{2}}(\boldsymbol{\mu}_{2.1},\boldsymbol{\Sigma}_{2.1};g^{(p_{2})})$
where
$\boldsymbol{\mu}_{2.1}=\boldsymbol{\mu}_{2}+\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}(\textbf{x}-\boldsymbol{\mu}_{1})$,
$\boldsymbol{\Sigma}_{2.1}=\boldsymbol{\Sigma}_{22}-\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12}$,
$g^{(p_{2})}(t)=\exp\\{-\frac{1}{2}(t+\delta_{1}(\textbf{x}))^{\beta}\\}$ and
$\delta_{1}(\textbf{x})=(\textbf{x}-\boldsymbol{\mu}_{1})^{\top}\boldsymbol{\Sigma}_{11}^{-1}(\textbf{x}-\boldsymbol{\mu}_{1})$
is the squared Mahalanobis distance (Gómez et al., 1998). Since sampling
directly from the marginal distribution of $\textbf{X}_{1}$ could be really
complicated, we will use a different approach that consists of drawing points
from the whole random vector of length $p$ and then approximate the moments
using Monte Carlo integration.

###  4.2 R function and Examples

The Algorithm 2 for the distributions mentioned in subsection 4.1 has been
implemented in the R function mvtelliptical, whose signature together with
default values is the following.

⬇

mvtelliptical(lower, upper=rep(Inf,length(lower)), mu=rep(0,length(lower)),

Sigma=diag(length(lower)), dist="Normal", nu=NULL, n=1e4, burn.in=0,

thinning=3)

The arguments lower and upper are the lower and upper truncation points of
length $p$, respectively, mu is the location parameter of length $p$, Sigma is
the $p\times p$ positive-definite scale matrix, nu is the additional parameter
or vector of parameters depending on the dgf $g$. The argument dist indicates
the distribution to be used. The parameters n, burn.in, and thinning are
related to the Monte Carlo approximation, where n is the number of samples to
be generated, burn.in is the number of samples to be discarded as burn-in
phase, and thinning is a factor for reducing autocorrelation between
observations.

Example 1

We illustrate how the method works considering a random vector of length 4
with truncated Student-$t$ distribution. In this example, the second variable
is not truncated, and the others are doubly truncated. The objective is to
study the performance of the estimates for the mean and the variance-
covariance elements obtained through Algorithm 2, considering a different
number of samples and thinning. After that, we compare those results with the
estimates from the R functions meanvarTMD available in package MomTrunc and
TT.moment from package TTmoment.

Figure 3 displays the boxplot for each element of the mean vector based on 100
estimates obtained through our proposal considering $n=10^{4}$ with thinning=1
and =3, $n=10^{5}$ with thinning=3, $n=3\times 10^{5}$ with thinning=1, and
$10^{6}$ samples with no thinning (=1). Also displays the results came from
the function meanvarTMD and the function TT.moment. The red dashed line
represents the median of the estimates achieved from the TT.moment function.
It is possible to observe that for the case of $n=10^{4}$, the estimates
obtained with no thinning have more variability than those with thinning=3
(observations with lower autocorrelation). The median of TT.moment estimates
is closer to the median of our method in most cases, except for $n=10^{4}$
with thinning=1. As expected, the variability in the estimates was reduced
when the sample size was increased. The distribution of the estimates from our
algorithm with $10^{5}$ samples and thinning=3 was similar to the distribution
considering $n=3\times 10^{5}$ and no thinning. Recall that both methods
needed to generate the same number of samples; the only difference here is
that the first one (thinning=3) will need less memory space than the other
one. The best results were obtained throughout TT.moment and meanvarTMD
functions. Those results are comparable with the estimates achieved from our
proposal with $n=10^{6}$ and no thinning.

Figure 3: Boxplot based on 100 estimates of the truncated mean. The red dashed
line represents the median of the estimates obtained from function TT.moment.
Figure 4: Boxplot based on 100 estimates of the variance-covariance elements.
The red dashed line represents the median of the estimates obtained from
function TT.moment.

Figure 4 shows the boxplot for the variance-covariance elements of the
truncated random vector considering each method described above. We noticed a
slight reduction in the variability of the estimates when considering a
thinning equal to 3. Another interesting fact is observed when we set
$n=10^{5}$ and thinning=3; in this case, it returned similar results than
estimate the covariances from MC with $3\times 10^{5}$ samples and no
thinning. The estimates achieved through our proposal considering $n=10^{6}$
are comparable with the results from TT.moment. The estimates obtained from
meanvarTMD are similar to those from MC with $n=10^{5}$ and thinning=3 in most
cases, except for $\sigma_{11},\sigma_{33}$, and $\sigma_{13}$. For these
parameters, our method showed better performance.

Example 2

In the previous example, it was observed that the estimates obtained from
Algorithm 2 with $n=10^{5}$ and thinning=3 are good enough to estimate the
mean and variance of a multivariate ($p=4$) variable with TMVT distribution,
even though the best results were gotten through the TT.moment function. In
this example, our goal is to analyze the execution time required for our
method and the functions meanvarTMD and TT.moment to estimate the first two
moments and the variance-covariance matrix of a $p$-variate random vector with
TMVT distribution considering $p=50,100,150$. In each case, we set 10%, 20%,
and 40% of the variables doubly truncated. The methods were run in a Windows
10 machine using R 4.0.3 on an Intel Core i7-7700 Processor with 3.60 GHz, and
32 GB of RAM.

Table 1 displays the median of the running time (in seconds) required for our
algorithm and functions meanvarTMD and TT.moment. For our proposal were
considered three scenarios $n=10^{4}$ with no thinning, $n=10^{4}$ with
thinning=3, and $n=10^{5}$ with thinning=3. The results are based on 100
simulations, and they were computed through the R function microbenchmark.
This table also shows the relative time computed, taking the time used by our
method with $n=10^{5}$ and thinning=3 as reference. We are going to refer to
this configuration as the “reference method". For our algorithm, we observed
that the time required to estimate the moments depends only on the number of
random observations sampled. Note that estimating the moments with $n=10^{4}$
took 3.50% of the time required for the reference method, and it is worth
mention that the number of samples needed for the first method is 3.33% the
number of samples used for the reference one. Our proposal with $n=10^{4}$ and
thinning=3 already needed 10% of the execution time used by the reference
method. Observe that the only case where meanvarTMD was faster than the
reference one is when the number of doubly truncated variables is equal to 5.
It also seems that the time needed by the meanvarTMD function depends only on
the number of doubly truncated variables. In all scenarios, the TT.moment
function is much more time-consuming, e.g., for a random vector of length
$p=100$ and 40 doubly truncated variables, it needed 28 times longer than the
reference method. An additional example can be found in Appendix A.2.

Table 1: Median of the execution time (in seconds) based on 100 simulations. Method | Measure | $p=50$ | $p=100$ | $p=150$
---|---|---|---|---
10% | 20% | 40% | 10% | 20% | 40% | 10% | 20% | 40%
$n=10^{4}$ | Median | 0.011 | 0.030 | 0.139 | 0.030 | 0.140 | 0.952 | 0.071 | 0.382 | 2.842
thinning = 1 | R.Time | 0.035 | 0.035 | 0.034 | 0.036 | 0.035 | 0.034 | 0.036 | 0.034 | 0.034
$n=10^{4}$ | Median | 0.031 | 0.084 | 0.404 | 0.085 | 0.405 | 2.820 | 0.199 | 1.118 | 8.461
thinning = 3 | R.Time | 0.100 | 0.100 | 0.100 | 0.100 | 0.100 | 0.100 | 0.101 | 0.100 | 0.100
$n=10^{5}$ | Median | 0.314 | 0.844 | 4.042 | 0.846 | 4.044 | 28.217 | 1.974 | 11.182 | 84.619
thinning = 3 | R.Time | - | - | - | - | - | - | - | - | -
meanvarTMD | Median | 0.118 | 4.102 | 49.189 | 3.781 | 48.681 | 367.243 | 21.209 | 157.179 | 1215.630
R.Time | 0.375 | 4.861 | 12.170 | 4.467 | 12.037 | 13.015 | 10.746 | 14.056 | 14.366
TT.moment | Median | 7.452 | 24.027 | 94.408 | 62.026 | 202.704 | 789.641 | 242.701 | 800.360 | 3081.367
R.Time | 23.767 | 28.0473 | 23.358 | 73.279 | 50.122 | 27.984 | 122.974 | 71.574 | 36.414

##  5 Application on Spatial Model for Censored Data

In this application we will consider the Gaussian spatial censored linear
(SCL) model defined by Lachos et al. (2017) and Ordoñez et al. (2018). In this
model, the data is generated from
$\textbf{Z}=\textbf{X}\boldsymbol{\beta}+\boldsymbol{\xi}$, with
$\boldsymbol{\xi}\sim\mathcal{N}_{p}(\textbf{0},\boldsymbol{\Sigma})$ and
$\boldsymbol{\Sigma}=[\mathrm{Cov}(\textbf{s}_{i},\textbf{s}_{j})]=\sigma^{2}\textbf{R}(\phi)+\tau^{2}\textbf{I}_{p}$.
It also has the particularity that the response variable Z is not fully
observed. Instead, it is observed $V_{i}$ and $C_{i}$ at each location, for
$i=1,\ldots,p$, where $C_{i}=0$ and $V_{i}=Z_{i}$ for an uncensored
observation $Z_{i}$, and if $C_{i}=1$ and $V_{i}=[V_{1i},V_{2i}]$ if $Z_{i}$
is censored or missing. Because of the difficulties in working directly with
the observed likelihood function, Lachos et al. (2017) suggested using an EM-
type algorithm to obtain the ML estimates of $\boldsymbol{\theta}$ considering
a parameterization $\boldsymbol{\Sigma}=\sigma^{2}\boldsymbol{\Psi}$, with
$\boldsymbol{\Psi}=\textbf{R}(\phi)+\nu^{2}\textbf{I}_{n}$ and
$\nu^{2}=\tau^{2}/\sigma^{2}$, to help with the identifiability of the
parameters. See also Diggle and Ribeiro (2007). Hence, the EM algorithm works
as it follows:

  * •

E-step: Let $\widehat{\boldsymbol{\theta}}^{(k)}$ be the current estimate of
$\boldsymbol{\theta}$, then the conditional expectation of the complete-data
log-likelihood without the constant is

$Q_{k}(\boldsymbol{\theta})=\mathbb{E}\left(\ell(\boldsymbol{\theta}|\textbf{Z}_{c})|\textbf{V},\textbf{C},\widehat{\boldsymbol{\theta}}^{(k)}\right)=-\frac{1}{2}\left[\log|\boldsymbol{\Psi}|+n\log\sigma^{2}+\frac{1}{\sigma^{2}}\widehat{A}^{(k)}\right],$

where
$\widehat{A}^{(k)}=\mathrm{tr}(\widehat{\textbf{Z}\textbf{Z}^{\top}}^{(k)}\boldsymbol{\Psi}^{-1})-2\widehat{\textbf{Z}}^{(k)\top}\boldsymbol{\Psi}^{-1}\textbf{X}\boldsymbol{\beta}+\boldsymbol{\beta}^{\top}\textbf{X}^{\top}\boldsymbol{\Psi}^{-1}\textbf{X}\boldsymbol{\beta}$.
Therefore, the E-step reduces only to the computation of
$\widehat{\textbf{Z}\textbf{Z}^{\top}}^{(k)}=\mathbb{E}\left(\textbf{Z}\textbf{Z}^{\top}|\textbf{V},\textbf{C},\widehat{\boldsymbol{\theta}}^{(k)}\right)$
and
$\widehat{\textbf{Z}}^{(k)}=\mathbb{E}\left(\textbf{Z}|\textbf{V},\textbf{C},\widehat{\boldsymbol{\theta}}^{(k)}\right)$.
In the traditional EM algorithm, we should now evaluate the conditional
expectations, which is possible through the R packages tmvtnorm or MomTrunc,
but it is computationally expensive when the proportion of censored
observations is non-negligible. An alternative is to use the MCEM algorithm,
which approximates the conditional expectations by using MC integration. For
the SCL model, the MCE-step is performed by estimating
$\widehat{\textbf{Z}\textbf{Z}}^{\top}$ and $\widehat{\textbf{Z}}$ through
Algorithm 2.

  * •

M-step: The conditional maximization step is carried out, and
$\widehat{\boldsymbol{\theta}}^{(k)}$ is updated by maximizing
$\widehat{Q}_{k}(\boldsymbol{\theta})$ over $\boldsymbol{\theta}$ to obtain a
new estimate $\widehat{\boldsymbol{\theta}}^{(k+1)}$, which leads to the
expressions:

$\displaystyle\widehat{\boldsymbol{\beta}}^{(k+1)}$ $\displaystyle=$
$\displaystyle\left(\textbf{X}^{\top}\widehat{\boldsymbol{\Psi}}^{-1(k)}\textbf{X}\right)^{-1}\textbf{X}^{\top}\widehat{\boldsymbol{\Psi}}^{-1(k)}\widehat{\textbf{Z}}^{(k)},$
$\displaystyle\widehat{\sigma}^{2(k+1)}$ $\displaystyle=$
$\displaystyle\frac{1}{n}\left[\mathrm{tr}\left(\widehat{\textbf{Z}\textbf{Z}^{\top}}^{(k)}\widehat{\boldsymbol{\Psi}}^{-1(k)}\right)-2\widehat{\textbf{Z}^{\top}}^{(k)}\widehat{\boldsymbol{\Psi}}^{-1(k)}\textbf{X}\widehat{\boldsymbol{\beta}}^{(k+1)}+\widehat{\boldsymbol{\beta}}^{\top(k+1)}\textbf{X}^{\top}\widehat{\boldsymbol{\Psi}}^{-1(k)}\textbf{X}\widehat{\boldsymbol{\beta}}^{(k+1)}\right],$
$\displaystyle\widehat{\boldsymbol{\alpha}}^{(k+1)}$ $\displaystyle=$
$\displaystyle\underset{\boldsymbol{\alpha}\in\mathbb{R}^{+}\times\mathbb{R}^{+}}{\operatorname{argmax}}\left(-\frac{1}{2}\log|\boldsymbol{\Psi}|-\frac{1}{2\widehat{\sigma}^{2(k+1)}}\left[\mathrm{tr}\left(\widehat{\textbf{Z}\textbf{Z}^{\top}}^{(k)}\boldsymbol{\Psi}^{-1}\right)-2\widehat{\textbf{Z}^{\top}}^{(k)}\boldsymbol{\Psi}^{-1}\textbf{X}\widehat{\boldsymbol{\beta}}^{(k+1)}\right.\right.$
$\displaystyle\left.\left.+\widehat{\boldsymbol{\beta}}^{\top(k+1)}\textbf{X}^{\top}\boldsymbol{\Psi}^{-1}\textbf{X}\widehat{\boldsymbol{\beta}}^{(k+1)}\right]\right),$

with $\boldsymbol{\alpha}=(\phi,\nu^{2})^{\top}$. Note that
$\widehat{\tau}^{2}$ can be recovered by
$\widehat{\tau}^{2(k+1)}=\widehat{\nu}^{2(k+1)}\widehat{\sigma}^{2(k+1)}$. An
efficient M-step can be easily accomplished by using, for instance, the roptim
package (Pan and Pan, 2020). In general, the estimates of
$\boldsymbol{\theta}$ may vary slightly around the maximum, with a variability
depending on the sample size used in the approximation. Hence, one possible
final estimate of the parameters may be computed as the mean of the estimates
after applying a burn-in and a thinning process.

###  5.1 Missouri Dioxin Contamination Data

The proposed MCEM algorithm will be applied to analyze the Missouri dioxin
contamination dataset available in CensSpatial package. The dataset contains
127 observations distributed in an area of $3600\times 65$ $m^{2}$ on the
shoulders of a country road located in Missouri, U.S.A. The observations
correspond to the level of contamination by dioxin
(2,3,7,8-tetrachlorodibenzo-p-dioxin or TCDD) at sampled points along the
road, where 43% of the observations (55 sites) were censored, falling below
some limit of detection, which ranges from 0.10 to 0.79 mg/kg. The spatial
directions are the $x$-direction (measured in 1/100 $ft$) and the
$y$-direction (in $ft$). Please, refer to Fridley and Dixon (2007) for more
details.

This dataset was firstly analyzed by Zirschky and Harris (1986), who concluded
that data appeared to be log-normally distributed. Hence, we fit the model
$\log(Z_{i})=\beta_{0}+\xi_{i}$, $i=1,\ldots,127$. The model parameters were
estimated using the MCEM algorithm and compared with the estimates from the
SAEM and EM algorithm. All methods were performed using 500 iterations and an
exponential correlation function to take into account the variation between
spatial points. For the MCEM algorithm, we evaluated four cases; in one of
those scenarios, it was considered linearly increasing sample sizes between
100 and 1000. Other scenarios considered constant sample sizes of 20, 5000,
and $10^{5}$. In order to use the SAEM algorithm, we set two configurations;
one draws points using the rmvtnorm function (from package tmvtnorm), and the
optimization procedure via optimx function (Nash et al., 2020). This method is
available in the CensSpatial package, and from now on, we refer to this
algorithm by SAEM. The second one draws points using the proposed slice
sampler, while the R function roptim executes the optimization procedure. We
will refer to the latter as SAEM-SS. Lastly, moments were computed using the
MomTrunc package for the EM algorithm. The functions used to estimate the
parameters via MCEM, SAEM-SS, and EM are available in the RcppCensSpatial
package.

Table 2: Missouri data - ML estimates and information criteria (AIC and BIC) obtained through MCEM, SAEM and EM algorithms considering the exponential correlation function. Algorithm | $n$ | $c$ | $\beta_{0}$ | $\sigma^{2}$ | $\phi$ | $\tau^{2}$ | Log-likelihood | AIC | BIC | Time (min)
---|---|---|---|---|---|---|---|---|---|---
MCEM | 20 |  | -2.355 | 6.577 | 14.702 | 0.213 | -143.128 | 294.257 | 305.633 | 0.936
$10^{2}$ \- $10^{3}$ | - | -2.402 | 6.808 | 15.095 | 0.207 | -143.108 | 294.216 | 305.592 | 1.819
$5000$ | - | -2.410 | 6.847 | 15.076 | 0.206 | -143.095 | 294.191 | 305.568 | 10.206
$10^{5}$ | - | -2.408 | 6.845 | 15.053 | 0.205 | -143.136 | 294.272 | 305.649 | 185.341
SAEM-SS | 20 | 0.25 | -2.332 | 6.312 | 15.109 | 0.214 | -143.153 | 294.307 | 305.683 | 0.676
SAEM | 20 | 0.25 | -2.014 | 4.858 | 14.206 | 0.245 | -143.840 | 295.681 | 307.057 | 6.079
$10^{5}$ | 1.00 | -2.010 | 4.829 | 14.136 | 0.245 | -143.865 | 295.729 | 307.106 | 9151.149
EM | - | - | -2.417 | 6.888 | 15.092 | 0.206 | -143.122 | 294.244 | 305.620 | 1661.472

The results of the ML estimates are shown in Table 2, where $n$ is the number
of samples considered to approximate the conditional mean, and $c$ indicates
the percentage of iterations without memory in the SAEM algorithm (Lachos et
al., 2017; Ordoñez et al., 2018). Final estimates for the MCEM and EM methods
were computed as the mean of the estimates at each iteration after applying a
burn-in of 250 and thinning of 3 observations, while the SAEM and SAEM-SS
estimates we only considered the estimates at the last iteration. We see that
estimates obtained from SAEM-SS are similar to MCEM estimates for $n=20$,
while the estimates for the EM algorithm are similar to MCEM estimates with
$n=5000$. The estimates obtained through MCEM and EM for the regression
coefficient $\beta_{0}$ were -2.400, while the SAEM algorithm estimated this
parameter equal to -2.010. Regarding to the spatial scaling parameter $\phi$,
it was around 15.05 and 14.10 for the MCEM and SAEM algorithms, respectively.
These values imply that for distances greater than 45 and 42 feet,
respectively, the correlation between two observations falls to less than
0.05. The estimates achieved from MCEM and SAEM methods for the partial sill
$\sigma^{2}$ and the nugget effect $\tau^{2}$ suggest that 97% and 95% of the
variability in data is explained by the spatial process, respectively. This
table also shows the maximized log-likelihood value, information criteria AIC
and BIC, and the running time in minutes. Based on information criteria, we
can conclude that MCEM with $n=5000$ best fits the Missouri dioxin
contamination data. Furthermore, it does not seem necessary to consider sample
sizes as large as $n=10^{5}$ because that configuration does not gain the
precision of the estimates and is more time-consuming.

Figure 5: Missouri data - Convergence of the parameter estimates via EM, MCEM,
and SAEM algorithm.

a. EM algorithm

b. MCEM algorithm

c. SAEM algorithm

Figure 5 shows the convergence graphs of the parameter estimates achieved from
MCEM, SAEM-SS, SAEM, and EM algorithms. Notice that the variability in the
estimates for MCEM decreases when the sample size increases from 100 to 1000
(aquamarine line). As expected, the estimates obtained from MCEM with $n=20$
(gray line) present more variability than the other three scenarios in which
we considered larger sample sizes, while MCEM with $n=10^{5}$ (black line)
reported the lowest variability in the estimates. The estimates of the
parameters computed through the EM algorithm present more variability than the
estimates from MCEM with $n=5000$ (red line), probably due to computational
stability of the numerical methods involved in the MomTrunc package; this is
why we decided to consider a burn-in and thinning procedure to compute the EM
final estimates.

##  6 Conclusions

This work describes an algorithm to generate random numbers from members of
the truncated elliptical family of distributions with a strictly decreasing
density generating function through a slice sampling algorithm and a Gibbs
sampler step. In addition, we presented an efficient approach to approximate
the first and the second moment for these distributions. We briefly introduce
the functions available in our R package relliptical in order to perform
sample generating and estimation of the first two moments. Simulation studies
were performed to investigate the properties of estimates and the robustness
of our algorithm. Moreover, we compared our approach with others available in
the R software (only for the normal and Student-$t$ case), where we showed
that our approach over-performed others in terms of precision and
computational time. We illustrate the usefulness of truncated moments on the
Missouri dioxin contamination dataset, where a spatial model for censored data
was fitted.

Future extensions of the work include the extension of this method to the
context of asymmetric multivariate elliptical distributions, so the fast
computation of their truncated moments may lead the way to proposed more
flexible and robust models relating censored models for mixed-effects models,
longitudinal data, spatial models, among others. Finally, results presented in
this paper can be reproduced through the R package relliptical, which is
available at CRAN for download.

## Acknowledgements

The research of Katherine A. L. Valeriano was supported by CAPES. Larissa A.
Matos acknowledges support from FAPESP-Brazil (Grant 2020/16713-0).

## References

  * Andersen et al. (2020) Andersen, M., R. Goedman, G. Grothendieck, S. Højsgaard, A. Pinkus, and G. Mazur (2020).  Ryacas: R interface to the YACAS Computer Algebra System.
  * Brent (2013) Brent, R. P. (2013).  Algorithms for minimization without derivatives.  Courier Corporation.
  * de Alencar et al. (2021) de Alencar, F. H., C. E. Galarza, L. A. Matos, and V. H. Lachos (2021).  Finite mixture modeling of censored and missing data using the multivariate skew-normal distribution.  Advances in Data Analysis and Classification, 1–37.
  * Delyon et al. (1999) Delyon, B., M. Lavielle, E. Moulines, et al. (1999).  Convergence of a stochastic approximation version of the EM algorithm.  The Annals of Statistics 27(1), 94–128.
  * Dempster et al. (1977) Dempster, A. P., N. M. Laird, and D. B. Rubin (1977).  Maximum likelihood from incomplete data via the EM algorithm.  Journal of the Royal Statistical Society: Series B (Methodological) 39(1), 1–22.
  * Diggle and Ribeiro (2007) Diggle, P. and P. Ribeiro (2007).  Model-based Geostatistics.  Springer.
  * Fang (2018) Fang, K. W. (2018).  Symmetric multivariate and related distributions.  CRC Press.
  * Fridley and Dixon (2007) Fridley, B. L. and P. Dixon (2007).  Data augmentation for a bayesian spatial model involving censored observations.  Environmetrics: The official journal of the International Environmetrics Society 18(2), 107–123.
  * Galarza et al. (2021) Galarza, C. E., R. Kan, and V. H. Lachos (2021).  Package ‘MomTrunc’.  R package version.
  * Galarza et al. (2021) Galarza, C. E., V. H. Lachos, and M. Bourguignon (2021).  A skew-t quantile regression for censored and missing data.  Stat 10(1), e379.
  * Galarza et al. (2021) Galarza, C. E., T.-I. Lin, W.-L. Wang, and V. H. Lachos (2021).  On moments of folded and truncated multivariate student-t distributions based on recurrence relations.  Metrika, 1–26.
  * Galarza et al. (2020) Galarza, C. E., L. A. Matos, and V. H. Lachos (2020).  Moments of the doubly truncated selection elliptical distributions with emphasis on the unified multivariate skew-$t$ distribution.  arXiv preprint arXiv:2007.14980.
  * Gelfand et al. (1992) Gelfand, A. E., A. F. Smith, and T.-M. Lee (1992).  Bayesian analysis of constrained parameter and truncated data problems using Gibbs sampling.  Journal of the American Statistical Association 87(418), 523–532.
  * Gómez et al. (1998) Gómez, E., M. Gomez-Viilegas, and J. M. Marín (1998).  A multivariate generalization of the power exponential family of distributions.  Communications in Statistics-Theory and Methods 27(3), 589–600.
  * Ho et al. (2015) Ho, H., T. Lin, W. Wang, A. Garay, V. Lachos, and M. Castro (2015).  R TTmoment package: sampling and calculating the first and second moments for the doubly truncated multivariate t distribution.
  * Ho et al. (2012) Ho, H. J., T.-I. Lin, H.-Y. Chen, and W.-L. Wang (2012).  Some results on the truncated multivariate t distribution.  Journal of Statistical Planning and Inference 142(1), 25–40.
  * Kan and Robotti (2017) Kan, R. and C. Robotti (2017).  On moments of folded and truncated multivariate normal distributions.  Journal of Computational and Graphical Statistics 26(4), 930–934.
  * Lachos et al. (2019) Lachos, V. H., L. A. Matos, L. M. Castro, and M.-H. Chen (2019).  Flexible longitudinal linear mixed models for multiple censored responses data.  Statistics in medicine 38(6), 1074–1102.
  * Lachos et al. (2017) Lachos, V. H., L. A. Matos, T. S. Barbosa, A. M. Garay, and D. K. Dey (2017).  Influence diagnostics in spatial models with censored response.  Environmetrics 28(7).
  * Matos et al. (2016) Matos, L. A., L. M. Castro, and V. H. Lachos (2016).  Censored mixed-effects models for irregularly observed repeated measures with applications to HIV viral loads.  Test 25(4), 627–653.
  * Matos et al. (2013) Matos, L. A., M. O. Prates, M.-H. Chen, and V. H. Lachos (2013).  Likelihood-based inference for mixed-effects models with censored response using the multivariate-t distribution.  Statistica Sinica 23(3), 1323–1345.
  * Morán-Vásquez and Ferrari (2019) Morán-Vásquez, R. A. and S. L. Ferrari (2019).  New results on truncated elliptical distributions.  Communications in Mathematics and Statistics, 1–15.
  * Muirhead (2009) Muirhead, R. J. (2009).  Aspects of multivariate statistical theory, Volume 197.  John Wiley & Sons.
  * Nash et al. (2020) Nash, J. C., R. Varadhan, and G. Grothendieck (2020).  Package ‘optimx’.
  * Neal (2003) Neal, R. M. (2003).  Slice sampling.  Annals of statistics, 705–741.
  * Ordoñez et al. (2018) Ordoñez, J. A., D. Bandyopadhyay, V. H. Lachos, and C. R. Cabral (2018).  Geostatistical estimation and prediction for censored responses.  Spatial statistics 23, 109–123.
  * Pan and Pan (2020) Pan, Y. and J. Pan (2020).  roptim: An R Package for General Purpose Optimization with C+.  R package version.
  * R Core Team (2021) R Core Team (2021).  R: A Language and Environment for Statistical Computing.  Vienna, Austria: R Foundation for Statistical Computing.
  * Robert and Casella (2010) Robert, C. P. and G. Casella (2010).  Introducing Monte Carlo Methods with R, Volume 18.  Springer.
  * Wei and Tanner (1990) Wei, G. C. and M. A. Tanner (1990).  A Monte Carlo implementation of the EM algorithm and the poor man’s data augmentation algorithms.  Journal of the American Statistical Association 85(411), 699–704.
  * Wilhelm (2015) Wilhelm, S. (2015).  Package ‘tmvtnorm’.  R journal.
  * Zirschky and Harris (1986) Zirschky, J. H. and D. J. Harris (1986).  Geostatistical analysis of hazardous waste site data.  Journal of Environmental Engineering 112(4), 770–784.

##  Appendix A Extra simulation results

This appendix contains additional information about the simulation results.

###  A.1 Sample autocorrelation

Figure 6: Sample autocorrelation plots of $X_{1}$ and $X_{2}$ sampled from the
bivariate truncated elliptical distributions in Figure 2.

###  A.2 Running Time to Compute Moments from Truncated Distributions

In this section, a complementary study of Example 2 (Section 4.1) was
conducted to examine the execution time required for our method in order to
estimate the first two moments and the variance-covariance matrix of a
$p$-variate random vector considering different distributions of the truncated
elliptical family, with $p=50$ and $100$. As in Example 2, for each case we
consider a 10%, 20%, and 40% of doubly truncated variables.

Table 3 shows the median of the execution time (in seconds) needed for
function mvtelliptical to compute the first two moments and the covariance
matrix. We considered a TMVN, a truncated contaminated normal with $\nu=1/2$
and $\rho=1/5$, a truncated Pearson VII with parameters $m=55$ and $\nu=3$, a
truncated slash with $\nu=2$ degrees of freedom, and a truncated power
exponential distribution with kurtosis $\beta=1/2$. For each case, our method
was applied setting $n=10^{4}$ and $10^{5}$ with a thinning$=3$. Notice that
the time needed by the algorithm for TMVN, TMVT, and truncated Pearson VII
distributions are similar and depend only on the number truncated variables
and samples used in the approximation. Our method requires more time to
compute moments from the truncated contaminated normal distribution when
compared to the latter results. This is because the algorithm uses a numerical
method to calculate the inverse of the dgf. Besides, it is interesting noting
that there is no time difference between computing the moments for a truncated
slash distribution with five or ten doubly truncated variables. This occurs
since the function used to approximate the integral on the dgf is more time-
consuming when $\nu+p/2-1$ is not an integer. Finally, the computation of the
moments for the truncated power exponential distribution required
approximately the same time for random vectors of equal length regardless of
the number of doubly truncated variables. For this case, the method samples
values for the whole vector, leading to no time difference.

Table 3: Median of the execution time (in seconds) based on 100 simulations. Distribution ($\nu$) | Sample size | $p=50$ | $p=100$
---|---|---|---
10% | 20% | 40% | 10% | 20% | 40%
Normal | $10^{4}$ | 0.028 | 0.083 | 0.399 | 0.081 | 0.400 | 2.888
$10^{5}$ | 0.285 | 0.840 | 3.999 | 0.805 | 4.003 | 28.892
Contaminated | $10^{4}$ | 0.071 | 0.118 | 0.440 | 0.120 | 0.442 | 2.928
Normal (1/2, 1/5) | $10^{5}$ | 0.706 | 1.180 | 4.405 | 1.192 | 4.415 | 29.286
Pearson VII (55, 3) | $10^{4}$ | 0.031 | 0.083 | 0.403 | 0.084 | 0.403 | 2.891
$10^{5}$ | 0.309 | 0.838 | 4.030 | 0.839 | 4.036 | 28.944
Slash (2) | $10^{4}$ | 0.202 | 0.202 | 0.548 | 0.200 | 0.549 | 3.113
$10^{5}$ | 2.020 | 2.026 | 5.481 | 1.997 | 5.489 | 31.160
Power | $10^{4}$ | 5.101 | 5.095 | 5.096 | 41.870 | 41.858 | 41.864
Exponential (1/2) | $10^{5}$ | 51.038 | 51.013 | 50.999 | 418.675 | 418.604 | 418.651

##  Appendix B The Multivariate Pearson VII Distribution

###  B.1 Marginal and conditional distributions

A random variable $\textbf{X}\in\mathbb{R}^{p}$ is said to have a multivariate
Pearson VII distribution with location parameter
$\boldsymbol{\mu}\in\mathbb{R}^{p}$, positive-definite scale matrix
$\boldsymbol{\Sigma}\in\mathbb{R}^{p\times p}$, extra parameters $m>p/2$ and
$\nu>0$, if its pdf is given by

$\displaystyle
f_{\textbf{X}}(\textbf{x})=\frac{\Gamma(m)}{(\pi\nu)^{p/2}\Gamma(m-p/2)}|\boldsymbol{\Sigma}|^{-1/2}\left(1+\frac{1}{\nu}(\textbf{x}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\textbf{x}-\boldsymbol{\mu})\right)^{-m},\quad\textbf{x}\in\mathbb{R}^{p}.$

The random vector X can also be represented as a scale mixture of normal (SMN)
distributions, i.e., $\textbf{X}=\boldsymbol{\mu}+U^{-1/2}\textbf{Z}$, where Z
has a $p$-variate normal distribution with mean $\textbf{0}\in\mathbb{R}^{p}$
and variance-covariance matrix $\boldsymbol{\Sigma}\in\mathbb{R}^{p\times p}$.
Here, $U$ follows Gamma distribution with scale parameter $m-p/2$ and rate
parameter $\nu/2$, where Z is independent of $U$. This implies that

$\textbf{X}|(U=u)\sim\mathcal{N}_{p}(\boldsymbol{\mu},u^{-1}\boldsymbol{\Sigma})\quad\mbox{and}\quad
U\sim\mathcal{G}(m-p/2,\nu/2).$

Therefore, the mean and the variance-covariance matrix of X are

$\displaystyle\mathbb{E}(\textbf{X})$ $\displaystyle=$
$\displaystyle\mathbb{E}(\mathbb{E}(\textbf{X}|U))=\mathbb{E}(\boldsymbol{\mu})=\boldsymbol{\mu},\quad
m>\frac{p+1}{2}.$ $\displaystyle\mathrm{\mathrm{C}ov}(\textbf{X})$
$\displaystyle=$
$\displaystyle\mathrm{Cov}(\mathbb{E}(\textbf{X}|U))+\mathbb{E}(\mathrm{Cov}(\textbf{X}|U))=\mathbb{E}(U^{-1})\boldsymbol{\Sigma}=\frac{\nu}{2m-p-2}\boldsymbol{\Sigma},\quad
m>\frac{p+2}{2}.\hskip 56.9055pt$

Now suppose that the vector X is partitioned into two random vectors
$\textbf{X}_{1}\in\mathbb{R}^{p_{1}}$ and
$\textbf{X}_{2}\in\mathbb{R}^{p_{2}}$, such that $p=p_{1}+p_{2}$, and consider
the partition of $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ used in
Proposition 4.1, i.e.,

$\displaystyle\textbf{X}=\left(\begin{array}[]{c}\textbf{X}_{1}\\\
\textbf{X}_{2}\end{array}\right),\quad\boldsymbol{\mu}=\left(\begin{array}[]{c}\boldsymbol{\mu}_{1}\\\
\boldsymbol{\mu}_{2}\end{array}\right)\quad\mbox{and}\quad\boldsymbol{\Sigma}=\left(\begin{array}[]{cc}\boldsymbol{\Sigma}_{11}&\boldsymbol{\Sigma}_{12}\\\
\boldsymbol{\Sigma}_{21}&\boldsymbol{\Sigma}_{22}\end{array}\right).$

First, notice that
$(\textbf{X}-\boldsymbol{\mu})^{\top}\boldsymbol{\Sigma}^{-1}(\textbf{X}-\boldsymbol{\mu})=\delta_{1}(\textbf{X}_{1})+\delta_{2.1}(\textbf{X}_{2.1})$,
where
$\delta_{1}(\textbf{X}_{1})=(\textbf{X}_{1}-\boldsymbol{\mu}_{1})^{\top}\boldsymbol{\Sigma}_{11}^{-1}(\textbf{X}_{1}-\boldsymbol{\mu}_{1})$,
$\delta_{2.1}(\textbf{X}_{2.1})=(\textbf{X}_{2}-\boldsymbol{\mu}_{2.1})^{\top}\boldsymbol{\Sigma}_{2.1}^{-1}(\textbf{X}_{2}-\boldsymbol{\mu}_{2.1})$,
$\boldsymbol{\mu}_{2.1}=\boldsymbol{\mu}_{2}+\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}(\textbf{X}_{1}-\boldsymbol{\mu}_{1})$
and
$\boldsymbol{\Sigma}_{2.1}=\boldsymbol{\Sigma}_{22}-\boldsymbol{\Sigma}_{21}\boldsymbol{\Sigma}_{11}^{-1}\boldsymbol{\Sigma}_{12}$.
By the results above, the marginal pdf of $\textbf{X}_{1}$ is given by

$\displaystyle f_{\textbf{X}_{1}}(\textbf{x}_{1})$ $\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{p_{2}}}f_{\textbf{X}}(\textbf{X})\mathrm{d}\textbf{x}_{2}=\frac{\Gamma(m)}{(\pi\nu)^{p/2}\Gamma(m-p/2)}|\boldsymbol{\Sigma}|^{-1/2}\int_{\mathbb{R}^{p_{2}}}\left(1+\frac{\delta_{1}(\textbf{x}_{1})}{\nu}+\frac{\delta_{2.1}(\textbf{x}_{2.1})}{\nu}\right)^{-m}\mathrm{d}\textbf{x}_{2}$
$\displaystyle=$
$\displaystyle\frac{\Gamma(m)}{(\pi\nu)^{p/2}\Gamma(m-p/2)}|\boldsymbol{\Sigma}|^{-1/2}\left(1+\frac{\delta_{1}(\textbf{x}_{1})}{\nu}\right)^{-m}\int_{\mathbb{R}^{p_{2}}}\left(1+\frac{\delta_{2.1}(\textbf{x}_{2.1})}{\nu+\delta_{1}(\textbf{x}_{1})}\right)^{-m}\mathrm{d}\textbf{x}_{2}$
$\displaystyle=$
$\displaystyle\frac{\Gamma(m-p_{2}/2)}{(\pi\nu)^{p_{1}/2}\Gamma(m-p/2)}|\boldsymbol{\Sigma}_{11}|^{-1/2}\left(1+\frac{\delta_{1}(\textbf{x}_{1})}{\nu}\right)^{-(m-p_{2}/2)},\quad\textbf{x}_{1}\in\mathbb{R}^{p_{1}}.$

Hence, the marginal distribution of $\textbf{X}_{1}$ is also Pearson VII
distributed with parameters $\boldsymbol{\mu}_{1}$,
$\boldsymbol{\Sigma}_{11}$, $m-p_{2}/2$ and $\nu$, i.e.,
$\textbf{X}_{1}\sim\mbox{PVII}_{p_{1}}(\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11},m-p_{2}/2,\nu)$.
On the other hand, the conditional pdf of
$\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x}_{1})$ is given by

$\displaystyle
f_{\textbf{X}_{2}|\textbf{X}_{1}}(\textbf{x}_{2}|\textbf{x}_{1})$
$\displaystyle=$
$\displaystyle\frac{f_{\textbf{X}}(\textbf{x}_{1},x_{2})}{f_{\textbf{X}_{1}}(\textbf{x}_{1})}$
$\displaystyle=$
$\displaystyle\frac{\Gamma(m)|\boldsymbol{\Sigma}_{2.1}|^{-1/2}}{(\pi(\nu+\delta_{1}(\textbf{x}_{1})))^{p_{2}/2}\Gamma(m-p_{2}/2)}\left(1+\frac{\delta_{2.1}(\textbf{x}_{2.1})}{\nu+\delta_{1}(\textbf{x}_{1})}\right)^{-m},\quad\textbf{x}_{1}\in\mathbb{R}^{p_{1}},\,\textbf{x}_{2}\in\mathbb{R}^{p_{2}}.$

Therefore, the conditional distribution has also a Pearson VII distribution
with parameters $\boldsymbol{\mu}_{2.1}$, $\boldsymbol{\Sigma}_{2.1}$, $m$ and
$\nu+\delta_{1}(\textbf{x}_{1})$, i.e.,
$\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x}_{1})\sim\mbox{PVII}_{p_{2}}(\boldsymbol{\mu}_{2.1},\boldsymbol{\Sigma}_{2.1},m,\nu+\delta_{1}(\textbf{x}_{1}))$.

###  B.2 Existence of its truncated moments

Let
$\textbf{X}\sim\mbox{PVII}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},m,\nu),m>p/2,\nu>0$,
and let ${A}\subseteq\mathbb{R}^{p}$ be a truncation region of interest. Then,
the expectation and the variance-covariance matrix of X given
$\textbf{X}\in{A}$ exist in the following cases:

  * •

If ${A}=\mathbb{R}^{p}$ or ${A}$ is unbounded, the vector is not truncated at
all, so the expectation exists for $m>(p+1)/2$ and the covariance matrix
exists for $m>(p+2)/2$, as usual.

  * •

If ${A}$ is bounded (all truncation points are finite), then
$\mathbb{E}(\textbf{X}|\textbf{X}\in{A})<\infty$ and
$\mathrm{Cov}(\textbf{X}|\textbf{X}\in{A})<\infty$ for all $m>p/2$, since the
distribution is bounded.

  * •

If X can be partitioned into two random variables
$\textbf{X}_{1}\in\mathbb{R}^{p_{1}}$ and
$\textbf{X}_{2}\in\mathbb{R}^{p_{2}}$ such that the truncation region
associated to $\textbf{X}_{1}$ (say, ${A}_{1}$) is bounded, from the last item
we have $\mathbb{E}(\textbf{X}_{1}|\textbf{X}\in{A})$ and
$\mathrm{Cov}(\textbf{X}_{1}|\textbf{X}\in{A})$ exist for all $m>p/2$ and
$\nu>0$. On the other hand, it follow from Fubini’s theorem that
$\mathbb{E}(\textbf{X}_{2}|\textbf{X}\in{A})$ will exist if and only if
$\mathbb{E}(\textbf{X}_{2}|\textbf{X}_{1})$ exists; this occurs for all
$m>(p_{2}+1)/2$. Note that $\mathbb{E}(\textbf{X}_{2}|\textbf{X}_{1})<\infty$
also implies that
$\mathrm{Cov}(\textbf{X}_{1},\textbf{X}_{2}|\textbf{X}\in{A})<\infty$.
Additionally, $\mathrm{Cov}(\textbf{X}_{2}|\textbf{X}\in{A})$ exists if and
only if $\mathrm{Cov}(\textbf{X}_{2}|\textbf{X}_{1})<\infty$, which holds for
all $m>(p_{2}+2)/2$.

###### Remark B.1

It is equivalent to say that $\mathbb{E}(\textbf{X}|\textbf{X}\in{A})$ exists
for all $m$, if at least one dimension containing a finite limit exists.
Besides, if at least two dimensions containing finite limits exist, we have
that $\mathrm{Cov}(\textbf{X}|\textbf{X}\in{A})$ exists for all $m>p/2$.

In order to illustrate the result, consider
$\textbf{X}\sim\mbox{PVII}_{2}(\boldsymbol{\mu},\boldsymbol{\Sigma},m,\nu)$,
with $\nu=1$, $\boldsymbol{\mu}=\textbf{0}$, and
$\boldsymbol{\Sigma}=\left(\begin{array}[]{cc}1&0.20\\\
0.20&1\end{array}\right)$. We are interesting to observe what happens with the
elements of $\mathbb{E}(\textbf{X}|\textbf{X}\in{A})$ and
$\mathrm{Cov}(\textbf{X}|\textbf{X}\in{A})$ for
${A}=\\{\textbf{x}\in\mathbb{R}^{2}:\textbf{a}<\textbf{x}<\textbf{b}\\}$ in
the following three scenarios:

  * a)

$m=2$, $\textbf{a}=(-0.80,-0.60)^{\top}$, $\textbf{b}=(\infty,\infty)^{\top}$;

  * b)

$m=1.40$, $\textbf{a}=(-0.80,-0.60)^{\top}$,
$\textbf{b}=(0.80,\infty)^{\top}$;

  * c)

$m=2$, $\textbf{a}=(-0.80,-0.60)^{\top}$, $\textbf{b}=(0.80,\infty)^{\top}$.

Figure 7 displays the trace evolution of the MC estimates for the mean and
variance-covariance elements $\mu_{1}$, $\mu_{2}$, $\sigma_{11}$,
$\sigma_{12}$ and $\sigma_{22}$ for each case. The red dashed line represents
the value for the parameter estimated via MC with $10^{6}$ samples, and we
refer to this value as the “real value".

For the first case, we have that $(p+1)/2=3/2<2=m$, then only the first moment
exists, i.e., $\mathbb{E}(\textbf{X}|\textbf{X}\in A)<\infty$. Therefore, we
observe in the first row of Figure 7 that only the estimates of $\mu_{1}$ and
$\mu_{2}$ converge to their real values as the sample size increase. In the
second scenario (middle row), we have that all elements converge except
$\sigma_{22}$. This happens because the truncation limits for the first
variable are finite and $m>(p_{2}+1)/2=1$. In the last case, scenario c),
convergence is attained for all parameters, since the condition
$m>(p_{2}+2)/2=3/2$ holds. Note that even with 2000 MC simulations there
exists a significant variability in the chains.

Figure 7: Trace plots of the evolution of the MC estimates for the mean and
variance-covariance elements of $\textbf{X}\mid(\textbf{X}\in A)$ under
scenarios a), b) and c). The red dashed line represents the true estimated
value computed using numerical methods

a. Two non-truncated variables, parameters $m=2$ and $\nu=1$.

b. One doubly truncated variable, parameters $m=1.40$ and $\nu=1$.

c. One doubly truncated variable, parameters $m=2$ and $\nu=1$.

##  Appendix C The Multivariate Slash Distribution

A random vector $\textbf{X}\in\mathbb{R}^{p}$ has multivariate slash
distribution with location parameter $\boldsymbol{\mu}\in\mathbb{R}^{p}$,
positive-definite scale matrix $\boldsymbol{\Sigma}\in\mathbb{R}^{p\times p}$
and $\nu>0$ degrees of freedom, denoted by
$\textbf{X}\sim\mbox{SL}_{p}(\boldsymbol{\mu},\boldsymbol{\Sigma},\nu)$, if
its pdf is given by

$f_{\textbf{X}}(\textbf{x})=\nu\int_{0}^{1}u^{\nu-1}\phi_{p}\left(\textbf{x};\boldsymbol{\mu},u^{-1}\boldsymbol{\Sigma}\right)\mathrm{d}u,\quad\textbf{x}\in\mathbb{R}^{p},$

where $\phi_{p}(\textbf{x};\boldsymbol{\mu},\boldsymbol{\Sigma})$ is the pdf
of a $p$-variate normal distribution with mean $\boldsymbol{\mu}$ and
covariance matrix $\boldsymbol{\Sigma}$. The pdf of a slash distribution can
be evaluated through numerical methods, e.g., using the R function integrate.
The random vector X can also be represented in the family of the SMN
distributions, this is, $\textbf{X}=\boldsymbol{\mu}+U^{-1/2}\textbf{Z}$,
where the random variables $U$ and Z are both independent and have
$\mathrm{Beta}(\nu,1)$ and $\mathcal{N}_{p}(\textbf{0},\boldsymbol{\Sigma})$
distributions, respectively. Therefore, the mean and variance-covariance
matrix of the random vector X are given by

$\displaystyle\mathbb{E}(\textbf{X})$ $\displaystyle=$
$\displaystyle\mathbb{E}\left(\mathbb{E}(\textbf{X}|U)\right)=\mathbb{E}(\boldsymbol{\mu})=\boldsymbol{\mu}.$
$\displaystyle\mathrm{Cov}(\textbf{X})$ $\displaystyle=$
$\displaystyle\mathrm{Cov}(\mathbb{E}(\textbf{X}|U))+\mathbb{E}(\mathrm{Cov}(\textbf{X}|U))=\mathbb{E}(U^{-1})\boldsymbol{\Sigma}=\frac{\nu}{\nu-1}\boldsymbol{\Sigma},\quad\nu>1.\hskip
85.35826pt$

Considering a partition in the same manner as used for the Pearson VII
distribution, the marginal pdf of $\textbf{X}_{1}$ is given by

$\displaystyle f_{\textbf{X}_{1}}(\textbf{x}_{1})$ $\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{p_{2}}}f_{\textbf{X}}(\textbf{x})d\textbf{x}_{2}=\int_{\mathbb{R}^{p_{2}}}\nu\int_{0}^{1}u^{\nu-1}\phi_{p}\left(\textbf{x};\boldsymbol{\mu},u^{-1}\boldsymbol{\Sigma}\right)\mathrm{d}u\,d\textbf{x}_{2}$
$\displaystyle=$
$\displaystyle\nu\int_{\mathbb{R}^{p_{2}}}\int_{0}^{1}u^{\nu-1}\phi_{p_{1}}\left(\textbf{x}_{1};\boldsymbol{\mu}_{1},u^{-1}\boldsymbol{\Sigma}_{11}\right)\phi_{p_{2}}\left(\textbf{x}_{2};\boldsymbol{\mu}_{2.1},u^{-1}\boldsymbol{\Sigma}_{2.1}\right)\mathrm{d}u\,\mathrm{d}\textbf{x}_{2}$
$\displaystyle\stackrel{{\scriptstyle\mbox{Fubini}}}{{=}}$
$\displaystyle\nu\int_{0}^{1}u^{\nu-1}\phi_{p_{1}}\left(\textbf{x}_{1};\boldsymbol{\mu}_{1},u^{-1}\boldsymbol{\Sigma}_{11}\right)\int_{\mathbb{R}^{p_{2}}}\phi_{p_{2}}\left(\textbf{x}_{2};\boldsymbol{\mu}_{2.1},u^{-1}\boldsymbol{\Sigma}_{2.1}\right)\mathrm{d}\textbf{x}_{2}\,\mathrm{d}u\hskip
85.35826pt$ $\displaystyle\Rightarrow
f_{\textbf{X}_{1}}(\textbf{x}_{1})=\nu\int_{0}^{1}u^{\nu-1}\phi_{p_{1}}\left(\textbf{x}_{1};\boldsymbol{\mu}_{1},u^{-1}\boldsymbol{\Sigma}_{11}\right)\mathrm{d}u.\hskip
227.62204pt$

Thus, $\textbf{X}_{1}\in\mathbb{R}^{p_{1}}$ follows a slash distribution with
location parameter $\boldsymbol{\mu}_{1}\in\mathbb{R}^{p_{1}}$, scale matrix
$\boldsymbol{\Sigma}_{11}\in\mathbb{R}^{p_{1}\times p_{1}}$ and $\nu>0$
degrees of freedom. On the other hand, the conditional pdf of
$\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x}_{1})$ is given by

$\displaystyle
f_{\textbf{X}_{2}|\textbf{X}_{1}}(\textbf{x}_{2}|\textbf{x}_{1})$
$\displaystyle=$
$\displaystyle\frac{f_{\textbf{X}}(\textbf{x}_{1},\textbf{x}_{2})}{f_{\textbf{X}_{1}}(\textbf{x}_{1})}=\frac{\nu}{f_{\textbf{X}_{1}}(\textbf{x}_{1})}\int_{0}^{1}u^{\nu-1}\phi_{p}\left(\textbf{x};\boldsymbol{\mu},u^{-1}\boldsymbol{\Sigma}\right)\mathrm{d}u$
$\displaystyle=$
$\displaystyle\frac{\nu}{f_{\textbf{X}_{1}}(\textbf{x}_{1})}\int_{0}^{1}u^{\nu-1}\phi_{p_{1}}\left(\textbf{x}_{1};\boldsymbol{\mu}_{1},u^{-1}\boldsymbol{\Sigma}_{11}\right)\phi_{p_{2}}\left(\textbf{x}_{2};\boldsymbol{\mu}_{2.1},u^{-1}\boldsymbol{\Sigma}_{2.1}\right)\mathrm{d}u.\hskip
56.9055pt$

Then, it is possible to notice that the Slash distribution is not closed under
conditioning. Furthermore, the pdf of
$\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x}_{1})$ belongs to the elliptical
family of distributions with dgf
$g(t)=\int_{0}^{1}u^{\nu+p/2-1}\exp\\{-u(t+\delta_{1}(\textbf{x}_{1}))/2\\}\mathrm{d}u$,
i.e.,
$\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x}_{1})\sim\mbox{E}\ell(\boldsymbol{\mu}_{2.1},\boldsymbol{\Sigma}_{2.1},\nu;g)$.
To determine the mean of the random vector
$\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x}_{1})$, we compute the conditional
expected value of the $i$th element of $\textbf{X}_{2}$ as follows

$\displaystyle\mathbb{E}(\textbf{X}_{2i}|\textbf{X}_{1}=\textbf{x}_{1})$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{p_{2}}}x_{2i}f_{\textbf{X}_{2}|\textbf{X}_{1}}(\textbf{x}_{2}|\textbf{x}_{1})\mathrm{d}\textbf{x}_{2}$
$\displaystyle=$
$\displaystyle\frac{\nu}{f_{\textbf{X}_{1}}(\textbf{x}_{1})}\int_{\mathbb{R}^{p_{2}}}x_{2i}\int_{0}^{1}u^{\nu-1}\phi_{p_{1}}\left(\textbf{x}_{1};\boldsymbol{\mu}_{1},u^{-1}\boldsymbol{\Sigma}_{11}\right)\phi_{p_{2}}\left(\textbf{x}_{2};\boldsymbol{\mu}_{2.1},u^{-1}\boldsymbol{\Sigma}_{2.1}\right)\mathrm{d}u\,\mathrm{d}\textbf{x}_{2}$
$\displaystyle\stackrel{{\scriptstyle\mbox{Fubini}}}{{=}}$
$\displaystyle\frac{\nu}{f_{\textbf{X}_{1}}(\textbf{x}_{1})}\int_{0}^{1}u^{\nu-1}\phi_{p_{1}}\left(\textbf{x}_{1};\boldsymbol{\mu}_{1},u^{-1}\boldsymbol{\Sigma}_{11}\right)\int_{\mathbb{R}^{p_{2}}}x_{2i}\phi_{p_{2}}\left(\textbf{x}_{2};\boldsymbol{\mu}_{2.1},u^{-1}\boldsymbol{\Sigma}_{2.1}\right)d\textbf{x}_{2}\,\mathrm{d}u$
$\displaystyle=$
$\displaystyle\frac{\mu_{2.1}^{(i)}\nu}{f_{\textbf{X}_{1}}(\textbf{x}_{1})}\int_{0}^{1}u^{\nu-1}\phi_{p_{1}}\left(\textbf{x}_{1};\boldsymbol{\mu}_{1},u^{-1}\boldsymbol{\Sigma}_{11}\right)\mathrm{d}u=\mu_{2.1}^{(i)},\quad\forall
i,\nu>0$

where $\mu_{2.1}^{(i)}$ represents the $i$th element of the vector
$\boldsymbol{\mu}_{2.1}$, and
$\mathbb{E}(\textbf{X}_{2}|\textbf{X}_{1}=\textbf{x}_{1})=\boldsymbol{\mu}_{2.1}$.
Now, to compute the elements of the variance-covariance matrix of the
conditional random vector, we first determine
$\mathbb{E}(X_{2i}X_{2j}|\textbf{X}_{1}=\textbf{x}_{1})$ for all
$i,j=1,\ldots,p_{2}$, as

$\displaystyle\mathbb{E}(X_{2i}X_{2j}|\textbf{X}_{1}=\textbf{x}_{1})$
$\displaystyle=$
$\displaystyle\int_{\mathbb{R}^{p_{2}}}x_{2i}x_{2j}f_{\textbf{X}_{2}|\textbf{X}_{1}}(\textbf{x}_{2}|\textbf{x}_{1})d\textbf{x}_{2}$
$\displaystyle=$
$\displaystyle\frac{\nu}{f_{\textbf{X}_{1}}(\textbf{x}_{1})}\int_{\mathbb{R}^{p_{2}}}x_{2i}x_{2j}\int_{0}^{1}u^{\nu-1}\phi_{p_{1}}\left(\textbf{x}_{1};\boldsymbol{\mu}_{1},u^{-1}\boldsymbol{\Sigma}_{11}\right)\phi_{p_{2}}\left(\textbf{x}_{2};\boldsymbol{\mu}_{2.1},u^{-1}\boldsymbol{\Sigma}_{2.1}\right)\mathrm{d}u\,\mathrm{d}\textbf{x}_{2}$
$\displaystyle\stackrel{{\scriptstyle\mbox{Fubini}}}{{=}}$
$\displaystyle\frac{\nu}{f_{\textbf{X}_{1}}(\textbf{x}_{1})}\int_{0}^{1}u^{\nu-1}\phi_{p_{1}}\left(\textbf{x}_{1};\boldsymbol{\mu}_{1},u^{-1}\boldsymbol{\Sigma}_{11}\right)\int_{\mathbb{R}^{p_{2}}}x_{2i}x_{2j}\phi_{p_{2}}\left(\textbf{x}_{2};\boldsymbol{\mu}_{2.1},u^{-1}\boldsymbol{\Sigma}_{2.1}\right)\mathrm{d}\textbf{x}_{2}\,\mathrm{d}u$
$\displaystyle=$
$\displaystyle\frac{\nu}{f_{\textbf{X}_{1}}(\textbf{x}_{1})}\int_{0}^{1}u^{\nu-1}\phi_{p_{1}}\left(\textbf{x}_{1};\boldsymbol{\mu}_{1},u^{-1}\boldsymbol{\Sigma}_{11}\right)\left(u^{-1}\sigma_{2.1}^{(ij)}+\mu_{2.1}^{(i)}\mu_{2.1}^{(j)}\right)\mathrm{d}u$
$\displaystyle=$
$\displaystyle\frac{\sigma_{2.1}^{(ij)}\nu}{f_{\textbf{X}_{1}}(\textbf{x}_{1})}\int_{0}^{1}u^{\nu-2}\phi_{p_{1}}\left(\textbf{x}_{1};\boldsymbol{\mu}_{1},u^{-1}\boldsymbol{\Sigma}_{11}\right)\mathrm{d}u+\mu_{2.1}^{(i)}\mu_{2.1}^{(j)}$
$\displaystyle=$
$\displaystyle\frac{\nu}{\nu-1}\left(\frac{\mbox{SL}_{p_{1}}(\textbf{x}_{1};\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11},\nu-1)}{\mbox{SL}_{p_{1}}(\textbf{x}_{1};\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11},\nu)}\right)\sigma_{2.1}^{(ij)}+\mu_{2.1}^{(i)}\mu_{2.1}^{(j)},\quad\nu>1,$

where $\sigma_{2.1}^{(ij)}$ is the $(i,j)$th element of the matrix
$\boldsymbol{\Sigma}_{2.1}$. From these results, we have that

$\displaystyle\mathrm{Cov}(X_{2i},X_{2j}|\textbf{X}_{1}=\textbf{x}_{1})=\frac{\nu}{\nu-1}\left(\frac{\mbox{SL}_{p_{1}}(\textbf{x}_{1};\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11},\nu-1)}{\mbox{SL}_{p_{1}}(\textbf{x}_{1};\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11},\nu)}\right)\sigma_{2.1}^{(ij)},\quad\nu>1.$

Therefore, the covariance matrix of the random vector
$\textbf{X}_{2}|(\textbf{X}_{1}=\textbf{x}_{1})$ will be given by

$\displaystyle\mathrm{Cov}(\textbf{X}_{2}|\textbf{X}_{1}=\textbf{x}_{1})=\frac{\nu}{\nu-1}\left(\frac{\mbox{SL}_{p_{1}}(\textbf{x}_{1};\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11},\nu-1)}{\mbox{SL}_{p_{1}}(\textbf{x}_{1};\boldsymbol{\mu}_{1},\boldsymbol{\Sigma}_{11},\nu)}\right)\boldsymbol{\Sigma}_{2.1},\quad\nu>1.$