`relliptical R`

packageThe `relliptical R`

package offers random numbers
generation from members of the truncated multivariate elliptical family
of distribution such as the truncated versions of the Normal, Student-t,
Pearson VII, Slash, Logistic, Kotz-type, among others. Particular
distributions can be provided by specifying the density generating
function. It also computes the first two moments (covariance matrix as
well) for some particular distributions. For more details see
(Valeriano, Galarza, and Matos 2021).

Next, we will show the functions available in the package.

The function `rtelliptical`

generates observations from a
truncated multivariate elliptical distribution with location parameter
`mu`

, scale matrix `Sigma`

, lower and upper
truncation points `lower`

and `upper`

via Slice
Sampling algorithm (Neal 2003) with Gibbs sampler (Robert and Casella
2010) steps. The `dist`

argument represents the truncated
distribution to be used. The values are `Normal`

,
`t`

, `PE`

, `PVII`

, `Slash`

,
and `CN`

, for the truncated Normal, Student-t, Power
Exponential, Pearson VII, Slash, and Contaminated Normal distributions,
respectively.

In the following example, we generate *n* = 10^{5}
samples from the truncated bivariate Normal distribution.

```
library(relliptical)
# Sampling from the Truncated Normal distribution
set.seed(1234)
= c(0, 1)
mu = matrix(c(3,0.6,0.6,3), 2, 2)
Sigma = c(-3, -3)
lower = c(3, 3)
upper = rtelliptical(n=1e5, mu, Sigma, lower, upper, dist="Normal")
sample1 head(sample1)
#> [,1] [,2]
#> [1,] 0.6643105 2.4005763
#> [2,] -1.3364441 -0.1756624
#> [3,] -0.1814043 1.7013605
#> [4,] -0.6841829 2.4750461
#> [5,] 2.0984490 0.1375868
#> [6,] -1.8796633 -1.2629126
library(ggplot2)
# Histogram and density for variable 1
= ggplot(data.frame(sample1), aes(x=X1)) +
f1 geom_histogram(aes(y=..density..), colour="black", fill="grey", bins=15) +
geom_density(colour="red") + labs(x=bquote(X[1]), y="Density")
# Histogram and density for variable 2
= ggplot(data.frame(sample1), aes(x=X2)) +
f2 geom_histogram(aes(y=..density..), colour="black", fill="grey", bins=15) +
geom_density(colour="red") + labs(x=bquote(X[2]), y="Density")
library(gridExtra)
grid.arrange(f1, f2, nrow=1)
```

This function also allows generating random numbers from other
truncated elliptical distributions not specified in the
`dist`

argument, by supplying the density generating function
(DGF) through arguments either `expr`

or `gFun`

.
The DGF must be a non-negative and strictly decreasing function on
(0,∞). 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`

and the
environment of `R`

. For instance, for the DGF
*g*(*t*) = *e*^{−t}, the user must
provide `expr = "exp(1)^(-t)"`

. See that the function must
depend only on variable *t*, and any additional parameter must be
passed as a fixed value. 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).

The following example draws random variates from a truncated
bivariate Logistic distribution, whose DGF is given by
*g*(*t*) = *e*^{−t}/(1+*e*^{−t})^{2}, *t* ≥ 0,
see (Fang, Kotz, and Ng 2018).

```
# Sampling from the Truncated Logistic distribution
= c(0, 0)
mu = matrix(c(1,0.70,0.70,1), 2, 2)
Sigma = c(-2, -2)
lower = c(3, 2)
upper # Sample autocorrelation with no thinning
set.seed(5678)
= rtelliptical(n=1e4, mu, Sigma, lower, upper, expr="exp(1)^(-t)/(1+exp(1)^(-t))^2")
sample2 tail(sample2)
#> [,1] [,2]
#> [9995,] -0.5639346 -1.4225548
#> [9996,] -0.4747796 -1.1890135
#> [9997,] -0.2596561 -0.1126482
#> [9998,] -1.1517285 -0.9875129
#> [9999,] -0.9268744 -1.5526844
#> [10000,] -1.1763390 -1.3098489
```

If it was no possible to generate random samples by passing a
character expression to `expr`

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

In the next example, we will draw samples from the truncated Kotz-type distribution, whose DGF is given by

*g*(*t*) = *t*^{N − 1}*e*^{−*r**ts}*, *t* ≥ 0, *r* > 0, *s* > 0, 2*N* + *p* > 2.

As required, this function is strictly decreasing when
(2−*p*)/2 < *N* ≤ 1, see (Fang, Kotz, and Ng 2018).

```
# Sampling from the Truncated Kotz-type distribution
set.seed(9876)
= c(0, 0)
mu = matrix(c(1,0.70,0.70,1), 2, 2)
Sigma = c(-2, -2)
lower = c(3, 2)
upper = rtelliptical(n=1e4, mu, Sigma, lower, upper, gFun=function(t){ t^(-1/2)*exp(-2*t^(1/4)) })
sample4 = ggplot(data.frame(sample4), aes(x=X1, y=X2)) + geom_point(size=0.50) +
f1 labs(x=expression(X[1]), y=expression(X[2]), subtitle="Kotz(2,1/4,1/2)")
library(ggExtra)
ggMarginal(f1, type="histogram", fill="grey")
```

Since the generating process uses an MCMC method, observations will be correlated, so it may be of interest to study some ACF plots. Now, we study the sample from the bivariate logistic distribution.

```
# Function for plotting the sample autocorrelation using ggplot2
= function(samples){
acf.plot = ncol(samples); n = nrow(samples); acf1 = list(p)
p for (i in 1:p){
= with(acf(samples[,i], plot=FALSE), data.frame(lag, acf))
bacfdf = ggplot(data=bacfdf, aes(x=lag,y=acf)) + geom_hline(aes(yintercept=0)) +
acf1[[i]] geom_segment(aes(xend=lag, yend=0)) + labs(x="Lag", y="ACF", subtitle=bquote(X[.(i)])) +
geom_hline(yintercept=c(qnorm(0.975)/sqrt(n),-qnorm(0.975)/sqrt(n)), colour="red", linetype="twodash")
}return (acf1)
}
grid.arrange(grobs=acf.plot(sample2), top="Sample ACF with no thinning", nrow=1)
```

Autocorrelation can be decimated by setting the `thinning`

argument. The thinning factor reduces the autocorrelation of the random
points in the Gibbs sampling process. Thinning consists in picking
separated points from the sample, at each k-th step. As natural, this
value must be an integer greater than or equal to 1.

```
# Sample autocorrelation with thinning = 3
set.seed(8768)
= rtelliptical(n=1e4, mu, Sigma, lower, upper, dist=NULL, expr="exp(1)^(-t)/(1+exp(1)^(-t))^2",
sample3 thinning=3)
grid.arrange(grobs=acf.plot(sample3), top="Sample ACF with thinning = 3", nrow=1)
```

For this purpose, we call the function `mvtelliptical()`

,
which returns the mean vector and variance-covariance matrix for some
specific truncated elliptical distributions. The argument
`dist`

sets the distribution to be used and accepts the same
values `Normal`

, `t`

, `PE`

,
`PVII`

, `Slash`

, and `CN`

as before.
Moments are computed through Monte Carlo method for the truncated
variables and using properties of the conditional expectation for the
non-truncated variables.

Next, we compute the moments for a random variable *X*
following a truncated 3-variate Student-t distribution with
*ν* = 0.8 degrees of freedom. We will consider two scenarios: a
first one with only one doubly truncated variable, and a second one with
two doubly truncated variables.

```
# Truncated Student-t distribution
set.seed(5678)
= c(0.1, 0.2, 0.3)
mu = matrix(data = c(1,0.2,0.3,0.2,1,0.4,0.3,0.4,1), nrow=length(mu), ncol=length(mu), byrow=TRUE)
Sigma # Example 1: considering nu = 0.80 and one doubly truncated variable
= c(-0.8, -Inf, -Inf)
a = c(0.5, 0.6, Inf)
b mvtelliptical(a, b, mu, Sigma, "t", 0.80)
#> $EY
#> [,1]
#> [1,] -0.11001805
#> [2,] -0.54278399
#> [3,] -0.01119847
#>
#> $EYY
#> [,1] [,2] [,3]
#> [1,] 0.13761136 0.09694152 0.04317817
#> [2,] 0.09694152 NaN NaN
#> [3,] 0.04317817 NaN NaN
#>
#> $VarY
#> [,1] [,2] [,3]
#> [1,] 0.12550739 0.03722548 0.04194614
#> [2,] 0.03722548 NaN NaN
#> [3,] 0.04194614 NaN NaN
# Example 2: considering nu = 0.80 and two doubly truncated variables
= c(-0.8, -0.70, -Inf)
a = c(0.5, 0.6, Inf)
b mvtelliptical(a, b, mu, Sigma, "t", 0.80) # By default n=1e4
#> $EY
#> [,1]
#> [1,] -0.08566441
#> [2,] 0.01563586
#> [3,] 0.19215627
#>
#> $EYY
#> [,1] [,2] [,3]
#> [1,] 0.126040187 0.005937196 0.01331868
#> [2,] 0.005937196 0.119761635 0.04700108
#> [3,] 0.013318682 0.047001083 1.14714388
#>
#> $VarY
#> [,1] [,2] [,3]
#> [1,] 0.118701796 0.007276632 0.02977964
#> [2,] 0.007276632 0.119517155 0.04399655
#> [3,] 0.029779636 0.043996554 1.11021985
```

As seen for the first scenario, some elements of the
variance-covariance matrix are shown as `NaN`

. Those are the
cases where the moment does not exist (yes, some elements of the
variance-covariance matrix may exist and others may not). It is well
know that for a Student-t distribution its second moment exist if
*ν* > 2, however, as studied by (Galarza et al. 2022), this
condition is relaxed as the number of dimensions containing only finite
truncation limits increases.

It is worth mention that the Student-t distribution with
*ν* > 0 degrees of freedom is a particular case of the Pearson
VII distribution with parameters *m* > *p*/2 and
*ν*^{*} > 0 when
*m* = (*ν*+*p*)/2 and
*ν*^{*} = *ν*.

Finally, for comparison purposes, we compute the moments for a
truncated Pearson VII distribution with parameters
*ν*^{*} = *ν* = 0.80 and
*m* = (*ν*+3)/2 = 1.90, which is equivalent to the
Student-t distribution above, and their moments must be nearly
equal.

```
# Truncated Pearson VII distribution
set.seed(9876)
= c(-0.8, -0.70, -Inf)
a = c(0.5, 0.6, Inf)
b = c(0.1, 0.2, 0.3)
mu = matrix(data = c(1,0.2,0.3,0.2,1,0.4,0.3,0.4,1), nrow=length(mu), ncol=length(mu), byrow=TRUE)
Sigma mvtelliptical(a, b, mu, Sigma, "PVII", c(1.90,0.80), n=1e6) # n=1e6 more precision
#> $EY
#> [,1]
#> [1,] -0.08558130
#> [2,] 0.01420611
#> [3,] 0.19166895
#>
#> $EYY
#> [,1] [,2] [,3]
#> [1,] 0.128348258 0.006903655 0.01420704
#> [2,] 0.006903655 0.121364742 0.04749544
#> [3,] 0.014207043 0.047495444 1.15156461
#>
#> $VarY
#> [,1] [,2] [,3]
#> [1,] 0.121024099 0.008119433 0.03061032
#> [2,] 0.008119433 0.121162929 0.04477257
#> [3,] 0.030610322 0.044772574 1.11482763
```

Fang, K. T., S. Kotz, and K. W. Ng. 2018. *Symmetric Multivariate and
Related Distributions*. Chapman; Hall/CRC.

Galarza, C. E., L. A. Matos, L. M. Castro, and V. H. Lachos. 2022.
“Moments of the Doubly Truncated Selection Elliptical Distributions with
Emphasis on the Unified Multivariate Skew-t Distribution.” *Journal
of Multivariate Analysis* 189: 104944. https://doi.org/10.1016/j.jmva.2021.104944.

Neal, R. M. 2003. “Slice Sampling.” *Annals of Statistics*,
705–41.

Robert, C. P., and G. Casella. 2010. *Introducing Monte Carlo Methods with R*. Vol. 18. New
York: Springer.

Valeriano, K. A., C. E. Galarza, and L. A. Matos. 2021. “Moments and
Random Number Generation for the Truncated Elliptical Family of
Distributions.” *arXiv Preprint arXiv:2112.09319*.