Type:	Package
Title:	Truncated Multivariate Normal and Student Distributions
Version:	2.3
Description:	A collection of functions to deal with the truncated univariate and multivariate normal and Student distributions, described in Botev (2017) <doi:10.1111/rssb.12162> and Botev and L'Ecuyer (2015) <doi:10.1109/WSC.2015.7408180>.
License:	GPL-3
BugReports:	https://github.com/lbelzile/TruncatedNormal/issues
Depends:	R (≥ 2.10)
Imports:	nleqslv, qrng, spacefillr, alabama, Rcpp (≥ 0.12.16)
LinkingTo:	Rcpp, RcppArmadillo
RoxygenNote:	7.3.2
VignetteBuilder:	knitr
Encoding:	UTF-8
Suggests:	knitr, rmarkdown, mvtnorm, carData, tinytest
NeedsCompilation:	yes
Packaged:	2024-07-08 15:35:59 UTC; lbelzile
Author:	Zdravko Botev [aut], Leo Belzile [aut, cre]
Maintainer:	Leo Belzile <belzilel@gmail.com>
Repository:	CRAN
Date/Publication:	2024-07-08 16:10:02 UTC

Truncated Normal Distribution Toolbox

Description

The routines include:

• generator of independent and identically distributed random vectors from the truncated univariate and multivariate distributions;

• (Quasi-) Monte Carlo estimator and a deterministic upper bound of the cumulative distribution function of the multivariate normal and Student distributions;

• algorithm for the accurate computation of the quantile function of the normal distribution in the extremes of its tails.

Author(s)

Leo Belzile and Z. I. Botev, email: botev@unsw.edu.au and web page: https://web.maths.unsw.edu.au/~zdravkobotev/

References

• Z. I. Botev (2017), The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting, Journal of the Royal Statistical Society, Series B, 79 (1), pp. 1–24.

• Z. I. Botev and P. L'Ecuyer (2015), Efficient Estimation and Simulation of the Truncated Multivariate Student-t Distribution, Proceedings of the 2015 Winter Simulation Conference, Huntington Beach, CA, USA

• Gibson G. J., Glasbey C. A., Elston D. A. (1994), Monte Carlo evaluation of multivariate normal integrals and sensitivity to variate ordering, In: Advances in Numerical Methods and Applications, pages 120–126

Cholesky matrix decomposition with GB ordering

Description

This function computes the Cholesky decomposition of a covariance matrix Sigma and returns a list containing the permuted bounds for integration. The prioritization of the variables follow the rule proposed in Gibson, Glasbey and Elston (1994) and reorder variables to have outermost variables with smallest expected values.

Usage

.cholpermGB(Sigma, l, u)

Arguments

Details

The list contains an integer vector perm with the indices of the permutation, which is such that Sigma(perm, perm) == L %*% t(L). The permutation scheme is described in Genz and Bretz (2009) in Section 4.1.3, p.37.

Value

a list with components

• L: Cholesky root

• l: permuted vector of lower bounds

• u: permuted vector of upper bounds

• perm: vector of integers with ordering of permutation

References

Genz, A. and Bretz, F. (2009). Computations of Multivariate Normal and t Probabilities, volume 105. Springer, Dordrecht.

Gibson G.J., Glasbey C.A. and D.A. Elton (1994). Monte Carlo evaluation of multivariate normal integrals and sensitivity to variate ordering. In: Dimon et al., Advances in Numerical Methods and Applications, WSP, pp. 120-126.

Cholesky matrix decomposition with GGE ordering

Description

This function computes the Cholesky decomposition of a covariance matrix Sigma and returns a list containing the permuted bounds for integration. The prioritization of the variables follow the rule proposed in Gibson, Glasbey and Elston (1994) and reorder variables to have outermost variables with smallest expected values.

Usage

.cholpermGGE(Sigma, l, u)

Arguments

Details

The list contains an integer vector perm with the indices of the permutation, which is such that Sigma(perm, perm) == L %*% t(L). The permutation scheme is described in Genz and Bretz (2009) in Section 4.1.3, p.37.

Value

a list with components

• L: Cholesky root

• l: permuted vector of lower bounds

• u: permuted vector of upper bounds

• perm: vector of integers with ordering of permutation

References

Genz, A. and Bretz, F. (2009). Computations of Multivariate Normal and t Probabilities, volume 105. Springer, Dordrecht.

Gibson G.J., Glasbey C.A. and D.A. Elton (1994). Monte Carlo evaluation of multivariate normal integrals and sensitivity to variate ordering. In: Dimon et al., Advances in Numerical Methods and Applications, WSP, pp. 120-126.

Multivariate normal density function

Description

This function returns the log-density for a multivariate Gaussian distribution. The data must be imputed as a matrix, using e.g., as.matrix, with each row representing an observation.

Usage

.dmvnorm_arma(x, mu, sigma, logd = FALSE)

Arguments

Value

density or log-density of the nrow(x) sample

Multivariate Student density function

Description

This function returns the log-density for a multivariate Student distribution. The data must be imputed as a matrix, using e.g., as.matrix, with each row representing an observation.

Usage

.dmvt_arma(x, mu, sigma, df, logd = FALSE)

Arguments

Value

density or log-density of the nrow(x) sample

Quantile of truncated Gaussian

Description

The function compute the pth quantile associated to the truncated standard Gaussian variate on the interval (l,u).

Usage

Phinv(p, l, u)

Arguments

Value

vector of quantiles

Cholesky decomposition for Gaussian distribution function with permutation

Description

This function computes the Cholesky decomposition of a covariance matrix Sigma and returns a list containing the permuted bounds for integration. The prioritization of the variables follows either the rule proposed in Gibson, Glasbey and Elston (1994), reorder variables to have outermost variables with smallest expected values. The alternative is the scheme proposed in Genz and Bretz (2009) that minimizes the variance of the truncated Normal variates.

Usage

cholperm(Sigma, l, u, method = c("GGE", "GB"))

Arguments

Details

The list contains an integer vector perm with the indices of the permutation, which is such that Sigma(perm, perm) == L %*% t(L). The permutation scheme is described in Genz and Bretz (2009) in Section 4.1.3, p.37.

Value

a list with components

• L: Cholesky root

• l: permuted vector of lower bounds

• u: permuted vector of upper bounds

• perm: vector of integers with ordering of permutation

References

Genz, A. and Bretz, F. (2009). Computations of Multivariate Normal and t Probabilities, volume 105. Springer, Dordrecht.

Gibson G.J., Glasbey C.A. and D.A. Elton (1994). Monte Carlo evaluation of multivariate normal integrals and sensitivity to variate ordering. In: Dimon et al., Advances in Numerical Methods and Applications, WSP, pp. 120-126.

Density function for the truncated multivariate normal distribution

Description

This function returns the (log)-density of a matrix x of observations lying in the interval [lb, ub].

Usage

dtmvnorm(
  x,
  mu,
  sigma,
  lb,
  ub,
  log = FALSE,
  type = c("mc", "qmc"),
  B = 10000,
  check = TRUE,
  ...
)

See Also

tmvnorm

Density function for the truncated multivariate Student distribution

Description

This function returns the (log)-density of a matrix x of observations lying in the interval [lb, ub].

Usage

dtmvt(
  x,
  mu,
  sigma,
  df,
  lb,
  ub,
  type = c("mc", "qmc"),
  log = FALSE,
  B = 10000,
  check = TRUE,
  ...
)

Arguments

See Also

tmvt

Calculate log of Gaussian distribution function accurately

Description

This function returns the probability of a standard Gaussian variate between the interval a and b, avoiding numerical overflow. The function is vectorized and is meant to be used only internally by the package TruncatedNormal.

Usage

lnNpr(a, b, check = TRUE)

Arguments

Value

a vector of log probability.

Latent membranous Lupus Nephritis dataset

Description

The data represents two clinical measurements (covariates), which are used to predict the occurrence of latent membranous lupus nephritis. The dataset consists of measurements on 55 patients of which 18 have been diagnosed with latent membranous lupus.

Format

a data frame with columns "response", "const", "x1" and "x2"

Source

The data were transcribed from Table 1, page 22, of Dyk and Meng (2001).

References

D. A. van Dyk and X.-L. Meng (2001) The art of data augmentation (with discussion). Journal of Computational and Graphical Statistics, volume 10, pages 1-50.

See Also

The dataset is used in the examples of mvrandn

Women wage dataset from Mroz (1987)

Description

The data are from the Panel Study of Income Dynamics (PSID) longitudinal study, 1976 wave. They give the number of work hours of married women along with socio-economic variables and the number of children.

Format

a data frame containing the following variables:

• whrs: hours of work

• kidslt6: number of children aged 5 and below years old in household

• kidsge6: number of children between age of 6 and 18 in household

• age: age (in years)

• educ: number of years in school

• hearn: hourly earnings

• exp: years of previous labor market experience

Source

W. Greene's website, accessed 17.12.2019 at <http://www.stern.nyu.edu/~wgreene/Text/Edition7/TableF5-1.csv>.

References

T. A. Mroz, 1987. The Sensitivity of an Empirical Model of Married Women's Hours of Work to Economic and Statistical Assumptions, Econometrica, 55(4), pp. 765-799

See Also

Mroz

Truncated multivariate normal cumulative distribution

Description

Computes an estimate and a deterministic upper bound of the probability Pr(l<X<u), where X is a zero-mean multivariate normal vector with covariance matrix \Sigma, that is, X is drawn from N(0,\Sigma). Infinite values for vectors u and l are accepted. The Monte Carlo method uses sample size n; the larger n, the smaller the relative error of the estimator.

Usage

mvNcdf(l, u, Sig, n = 1e+05)

Arguments

Details

Suppose you wish to estimate Pr(l<AX<u), where A is a full rank matrix and X is drawn from N(\mu,\Sigma), then you simply compute Pr(l-A\mu<AY<u-A\mu), where Y is drawn from N(0, A\Sigma A^\top).

Value

a list with components

• prob: estimated value of probability Pr(l<X<u)

• relErr: estimated relative error of estimator

• upbnd: theoretical upper bound on true Pr(l<X<u)

Note

For small dimensions, say d<50, better accuracy may be obtained by using the (usually slower) quasi-Monte Carlo version mvNqmc of this algorithm.

Author(s)

Zdravko I. Botev

References

Z. I. Botev (2017), The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting, Journal of the Royal Statistical Society, Series B, 79 (1), pp. 1–24.

See Also

mvNqmc, mvrandn

Examples

d <- 15; l <- 1:d; u <- rep(Inf, d);
Sig <- matrix(rnorm(d^2), d, d)*2; Sig=Sig %*% t(Sig)
mvNcdf(l, u, Sig, 1e4) # compute the probability 

Truncated multivariate normal cumulative distribution (quasi-Monte Carlo)

Description

Computes an estimate and a deterministic upper bound of the probability Pr(l<X<u), where X is a zero-mean multivariate normal vector with covariance matrix \Sigma, that is, X is drawn from N(0,\Sigma). Infinite values for vectors u and l are accepted. The Monte Carlo method uses sample size n: the larger n, the smaller the relative error of the estimator.

Usage

mvNqmc(l, u, Sig, n = 1e+05)

Arguments

Details

Suppose you wish to estimate Pr(l<AX<u), where A is a full rank matrix and X is drawn from N(\mu,\Sigma), then you simply compute Pr(l-A\mu<AY<u-A\mu), where Y is drawn from N(0, A\Sigma A^\top).

Value

a list with components

• prob: estimated value of probability Pr(l<X<u)

• relErr: estimated relative error of estimator

• upbnd: theoretical upper bound on true Pr(l<X<u)

Note

This version uses a Quasi Monte Carlo (QMC) pointset of size ceiling(n/12) and estimates the relative error using 12 independent randomized QMC estimators. QMC is slower than ordinary Monte Carlo, but is also likely to be more accurate when d<50. For high dimensions, say d>50, you may obtain the same accuracy using the (typically faster) mvNcdf.

Author(s)

Zdravko I. Botev

References

Z. I. Botev (2017), The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting, Journal of the Royal Statistical Society, Series B, 79 (1), pp. 1–24.

See Also

mvNcdf, mvrandn

Examples

d <- 15
l <- 1:d
u <- rep(Inf, d)
Sig <- matrix(rnorm(d^2), d, d)*2
Sig <- Sig %*% t(Sig)
mvNqmc(l, u, Sig, 1e4) # compute the probability

Truncated multivariate student cumulative distribution

Description

Computes an estimator of the probability Pr(l<X<u), where X is a centered multivariate student vector with scale matrix Sig and degrees of freedom df. Infinite values for vectors u and l are accepted.

Usage

mvTcdf(l, u, Sig, df, n = 1e+05)

Arguments

Details

Monte Carlo method uses sample size n; the larger the n, the smaller the relative error of the estimator;

Value

a list with components

• prob: estimated value of probability Pr(l<X<u)

• relErr: estimated relative error of estimator

• upbnd: theoretical upper bound on true Pr(l<X<u)

Note

If you want to estimate Pr(l<Y<u), where Y follows a Student distribution with df degrees of freedom, location vector m and scale matrix Sig, then use mvTqmc(Sig, l - m, u - m, nu, n).

Author(s)

Matlab code by Zdravko Botev, R port by Leo Belzile

References

Z. I. Botev (2017), The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting, Journal of the Royal Statistical Society, Series B, 79 (1), pp. 1–24

Z. I. Botev and P. L'Ecuyer (2015), Efficient probability estimation and simulation of the truncated multivariate Student-t distribution, Proceedings of the 2015 Winter Simulation Conference, pp. 380-391

See Also

mvTqmc, mvrandt, mvNqmc, mvrandn

Examples

 d <- 15; nu <- 30;
 l <- rep(2, d); u <- rep(Inf, d);
 Sig <- 0.5 * matrix(1, d, d) + 0.5 * diag(1, d);
 est <- mvTcdf(l, u, Sig, nu, n = 1e4)
 # mvtnorm::pmvt(lower = l, upper = u, df = nu, sigma = Sig)
## Not run: 
d <- 5
Sig <- solve(0.5*diag(d)+matrix(0.5, d,d))
# mvtnorm::pmvt(lower = rep(-1,d), upper = rep(Inf, d), df = 10, sigma = Sig)[1]
mvTcdf(rep(-1, d), u = rep(Inf, d), Sig = Sig, df = 10, n=1e4)$prob

## End(Not run)

Truncated multivariate student cumulative distribution (QMC version)

Description

Computes an estimator of the probability Pr(l<X<u), where X is a zero-mean multivariate student vector with scale matrix Sig and degrees of freedom df. Infinite values for vectors u and l are accepted.

Usage

mvTqmc(l, u, Sig, df, n = 1e+05)

Arguments

Details

This version uses a Quasi Monte Carlo (QMC) pointset of size ceiling(n/12) and estimates the relative error using 12 independent randomized QMC estimators; QMC is slower than ordinary Monte Carlo (see mvTcdf), but is also likely to be more accurate when d<50.

Value

a list with components

• prob: estimated value of probability Pr(l<X<u)

• relErr: estimated relative error of estimator

• upbnd: theoretical upper bound on true Pr(l<X<u)

Note

If you want to estimate Pr(l<Y<u), where Y follows a Student distribution with df degrees of freedom, location vector m and scale matrix Sig, then use mvTqmc(Sig, l - m, u - m, nu, n).

Author(s)

Matlab code by Zdravko I. Botev, R port by Leo Belzile

References

Z. I. Botev (2017), The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting, Journal of the Royal Statistical Society, Series B, 79 (1), pp. 1–24

Z. I. Botev and P. L'Ecuyer (2015), Efficient probability estimation and simulation of the truncated multivariate Student-t distribution, Proceedings of the 2015 Winter Simulation Conference, pp. 380-391

See Also

mvTcdf, mvrandt, mvNqmc, mvrandn

Examples

d <- 25; nu <- 30;
l <- rep(1, d) * 5; u <- rep(Inf, d);
Sig <- 0.5 * matrix(1, d, d) + 0.5 * diag(d);
est <- mvTqmc(l, u, Sig, nu, n = 1e4)
## Not run: 
d <- 5
Sig <- solve(0.5*diag(d)+matrix(0.5, d,d))
## mvtnorm::pmvt(lower = rep(-1,d), upper = rep(Inf, d), df = 10, sigma = Sig)[1]
mvTqmc(rep(-1, d), u = rep(Inf, d), Sig = Sig, df = 10, n=1e4)$prob

## End(Not run)

Truncated multivariate normal generator

Description

Simulate n independent and identically distributed random vectors from the d-dimensional N(0,\Sigma) distribution (zero-mean normal with covariance \Sigma) conditional on l<X<u. Infinite values for l and u are accepted.

Usage

mvrandn(l, u, Sig, n, mu = NULL)

Arguments

Details

• Bivariate normal: Suppose we wish to simulate a bivariate X from N(\mu,\Sigma), conditional on X_1-X_2<-6. We can recast this as the problem of simulation of Y from N(0,A\Sigma A^\top) (for an appropriate matrix A) conditional on l-A\mu < Y < u-A\mu and then setting X=\mu+A^{-1}Y. See the example code below.

• Exact posterior simulation for Probit regression: Consider the Bayesian Probit Regression model applied to the lupus dataset. Let the prior for the regression coefficients \beta be N(0,\nu^2 I). Then, to simulate from the Bayesian posterior exactly, we first simulate Z from N(0,\Sigma), where \Sigma=I+\nu^2 X X^\top, conditional on Z\ge 0. Then, we simulate the posterior regression coefficients, \beta, of the Probit regression by drawing (\beta|Z) from N(C X^\top Z,C), where C^{-1}=I/\nu^2+X^\top X. See the example code below.

Value

a d by n matrix storing the random vectors, X, drawn from N(0,\Sigma), conditional on l<X<u;

Note

The algorithm may not work or be very inefficient if \Sigma is close to being rank deficient.

See Also

mvNqmc, mvNcdf

Examples

 # Bivariate example.

 Sig <- matrix(c(1,0.9,0.9,1), 2, 2);
 mu <- c(-3,0); l <- c(-Inf,-Inf); u <- c(-6,Inf);
 A <- matrix(c(1,0,-1,1),2,2);
 n <- 1e3; # number of sampled vectors
 Y <- mvrandn(l - A %*% mu, u - A %*% mu, A %*% Sig %*% t(A), n);
 X <- rep(mu, n) + solve(A, diag(2)) %*% Y;
 # now apply the inverse map as explained above
 plot(X[1,], X[2,]) # provide a scatterplot of exactly simulated points
## Not run: 
# Exact Bayesian Posterior Simulation Example.

data("lupus"); # load lupus data
Y = lupus[,1]; # response data
X = lupus[,-1]  # construct design matrix
m=dim(X)[1]; d=dim(X)[2]; # dimensions of problem
 X=diag(2*Y-1) %*%X; # incorporate response into design matrix
 nu=sqrt(10000); # prior scale parameter
 C=solve(diag(d)/nu^2+t(X)%*%X);
 L=t(chol(t(C))); # lower Cholesky decomposition
 Sig=diag(m)+nu^2*X %*% t(X); # this is covariance of Z given beta
 l=rep(0,m);u=rep(Inf,m);
 est=mvNcdf(l,u,Sig,1e3);
 # estimate acceptance probability of Crude Monte Carlo
 print(est$upbnd/est$prob)
 # estimate the reciprocal of acceptance probability
 n=1e4 # number of iid variables
 z=mvrandn(l,u,Sig,n);
 # sample exactly from auxiliary distribution
 beta=L %*% matrix(rnorm(d*n),d,n)+C %*% t(X) %*% z;
 # simulate beta given Z and plot boxplots of marginals
 boxplot(t(beta))
 # plot the boxplots of the marginal
 # distribution of the coefficients in beta
 print(rowMeans(beta)) # output the posterior means
 
## End(Not run)

Random number generation from multivariate truncated Student distribution

Description

Random number generation from multivariate truncated Student distribution

Usage

mvrandt(l, u, Sig, df, n, mu = NULL)

Arguments

Value

a d by n matrix

Author(s)

Matlab code by Zdravko Botev, R port by Leo Belzile

References

Z. I. Botev and P. L'Ecuyer (2015), Efficient probability estimation and simulation of the truncated multivariate Student-t distribution, Proceedings of the 2015 Winter Simulation Conference, pp. 380-391,

Examples

## Not run: 
d <- 60L; n <- 1e3;
Sig <- 0.9 * matrix(1, d, d) + 0.1 * diag(d);
l <- (1:d)/d * 4; u <- l+2; df <- 10;
X <- mvrandt(l,u,Sig,df,n)
stopifnot(all(X>l))
stopifnot(all(X<u))

## End(Not run)

Normal quantile function (high precision)

Description

Computes with tail-precision the quantile function of the standard normal distribution at 0\le p\le 1, and truncated to the interval [l,u]. Infinite values for vectors l and u are accepted.

Usage

norminvp(p, l, u)

Arguments

Details

Suppose we wish to simulate a random variable Z drawn from N(\mu,\sigma^2) and conditional on l<Z<u using the inverse transform method. To achieve this, first compute X=norminvp(runif(1),(l-mu)/sig,(u-mu)/sig) and then set Z=mu+sig*X

Value

quantile value of the truncated normal distribution.

Note

If you wish to simulate truncated normal variables fast, use trandn. Using norminvp is advisable only when needed, for example, in quasi-Monte Carlo or antithetic sampling, where the inverse transform method is unavoidable.

Author(s)

Zdravko I. Botev

References

Z. I. Botev (2017), The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting, Journal of the Royal Statistical Society, Series B, 79 (1), pp. 1–24.

See Also

trandn

Examples

 d <- 150 # simulate via inverse transform method
 norminvp(runif(d),l = 1:d, u = rep(Inf, d))

Distribution function of the multivariate normal distribution for arbitrary limits

Description

This function computes the distribution function of a multivariate normal distribution vector for an arbitrary rectangular region [lb, ub]. pmvnorm computes an estimate and the value is returned along with a relative error and a deterministic upper bound of the distribution function of the multivariate normal distribution. Infinite values for vectors u and l are accepted. The Monte Carlo method uses sample size n: the larger the sample size, the smaller the relative error of the estimator.

Usage

pmvnorm(
  mu,
  sigma,
  lb = -Inf,
  ub = Inf,
  B = 10000,
  type = c("mc", "qmc"),
  check = TRUE,
  ...
)

Arguments

Author(s)

Zdravko I. Botev, Leo Belzile (wrappers)

References

Z. I. Botev (2017), The normal law under linear restrictions: simulation and estimation via minimax tilting, Journal of the Royal Statistical Society, Series B, 79 (1), pp. 1–24.

See Also

pmvnorm

Examples

#From mvtnorm
mean <- rep(0, 5)
lower <- rep(-1, 5)
upper <- rep(3, 5)
corr <- matrix(0.5, 5, 5) + diag(0.5, 5)
prob <- pmvnorm(lb = lower, ub = upper, mu = mean, sigma = corr)
stopifnot(pmvnorm(lb = -Inf, ub = 3, mu = 0, sigma = 1) == pnorm(3))

Distribution function of the multivariate Student distribution for arbitrary limits

Description

This function computes the distribution function of a multivariate normal distribution vector for an arbitrary rectangular region [lb, ub]. pmvt computes an estimate and the value is returned along with a relative error and a deterministic upper bound of the distribution function of the multivariate normal distribution. Infinite values for vectors u and l are accepted. The Monte Carlo method uses sample size n: the larger the sample size, the smaller the relative error of the estimator.

Usage

pmvt(
  mu,
  sigma,
  df,
  lb = -Inf,
  ub = Inf,
  type = c("mc", "qmc"),
  B = 10000,
  check = TRUE,
  ...
)

Arguments

Author(s)

Matlab code by Zdravko I. Botev, R port by Leo Belzile

References

Z. I. Botev and P. L'Ecuyer (2015), Efficient probability estimation and simulation of the truncated multivariate Student-t distribution, Proceedings of the 2015 Winter Simulation Conference, pp. 380-391

Examples

d <- 15; nu <- 30;
l <- rep(2, d); u <- rep(Inf, d);
sigma <- 0.5 * matrix(1, d, d) + 0.5 * diag(1, d);
est <- pmvt(lb = l, ub = u, sigma = sigma, df = nu)
# mvtnorm::pmvt(lower = l, upper = u, df = nu, sigma = sigma)
## Not run: 
d <- 5
sigma <- solve(0.5 * diag(d) + matrix(0.5, d, d))
# mvtnorm::pmvt(lower = rep(-1,d), upper = rep(Inf, d), df = 10, sigma = sigma)[1]
pmvt(lb = rep(-1, d), ub = rep(Inf, d), sigma = sigma, df = 10)

## End(Not run)

Cumulative distribution function of the truncated multivariate normal distribution.

Description

This function returns the (log)-distribution function of a matrix q of observations lying in the interval [lb, ub].

Usage

ptmvnorm(
  q,
  mu,
  sigma,
  lb,
  ub,
  log = FALSE,
  type = c("mc", "qmc"),
  B = 10000,
  check = TRUE,
  ...
)

See Also

tmvnorm

Cumulative distribution function of the truncated multivariate Student distribution.

Description

This function returns the (log)-distribution function of a matrix q of observations lying in the interval [lb, ub].

Usage

ptmvt(
  q,
  mu,
  sigma,
  df,
  lb,
  ub,
  type = c("mc", "qmc"),
  log = FALSE,
  check = TRUE,
  B = 10000,
  ...
)

Arguments

See Also

tmvt

Quantile function using the inversion method

Description

If p is a matrix, the arguments are recycled. If mu or sd are not specified they assume the default values of 0 and 1, respectively.

Usage

qtnorm(p, mu, sd, lb, ub)

See Also

tnorm

Random number generator for the truncated multivariate normal distribution.

Description

This function returns a matrix of draws from a multivariate normal distribution truncated on the interval [lb, ub].

Usage

rtmvnorm(n, mu, sigma, lb, ub, check = TRUE, ...)

See Also

tmvnorm

Random number generator for the truncated multivariate Student distribution.

Description

This function returns a matrix of draws from a multivariate Student distribution truncated on the interval [lb, ub].

Usage

rtmvt(n, mu, sigma, df, lb, ub, check = TRUE, ...)

Arguments

See Also

tmvt

Vectorized random number generation from the univariate truncated normal distribution

Description

If n is 1, the function returns a vector rather than a matrix. If mu or sd are not specified they assume the default values of 0 and 1, respectively.

Usage

rtnorm(n, mu, sd, lb, ub, method = c("fast", "invtransfo"))

See Also

tnorm

Multivariate truncated normal distribution

Description

Density, distribution function and random generation for the multivariate truncated normal distribution with mean vector mu, covariance matrix sigma, lower truncation limit lb and upper truncation limit ub. The truncation limits can include infinite values. The Monte Carlo (type = "mc") uses a sample of size B, while the quasi Monte Carlo (type = "qmc") uses a pointset of size ceiling(n/12) and estimates the relative error using 12 independent randomized QMC estimators.

Arguments

Value

dtmvnorm gives the density, ptmvnorm and pmvnorm give the distribution function of respectively the truncated and multivariate Gaussian distribution and rtmvnorm generate random deviates.

Usage

dtmvnorm(x, mu, sigma, lb, ub, type = c("mc", "qmc"), log = FALSE, B = 1e4)
ptmvnorm(q, mu, sigma, lb, ub, type = c("mc", "qmc"), log = FALSE, B = 1e4)
rtmvnorm(n, mu, sigma, lb, ub)

Author(s)

Zdravko I. Botev, Leo Belzile (wrappers)

References

Z. I. Botev (2017), The normal law under linear restrictions: simulation and estimation via minimax tilting, Journal of the Royal Statistical Society, Series B, 79 (1), pp. 1–24.

Examples

d <- 4; lb <- rep(0, d)
mu <- runif(d)
sigma <- matrix(0.5, d, d) + diag(0.5, d)
samp <- rtmvnorm(n = 10, mu = mu, sigma = sigma, lb = lb)
loglik <- dtmvnorm(samp, mu = mu, sigma = sigma, lb = lb, log = TRUE)
cdf <- ptmvnorm(samp, mu = mu, sigma = sigma, lb = lb, log = TRUE, type = "q")

# Exact Bayesian Posterior Simulation Example
# Vignette, example 5
## Not run: 
data("lupus"); # load lupus data
Y <- lupus[,1]; # response data
X <- as.matrix(lupus[,-1])  # construct design matrix
n <- nrow(X)
d <- ncol(X)
X <- diag(2*Y-1) %*% X; # incorporate response into design matrix
nusq <- 10000; # prior scale parameter
C <- solve(diag(d)/nusq + crossprod(X))
sigma <- diag(n) + nusq*tcrossprod(X) # this is covariance of Z given beta
est <- pmvnorm(sigma = sigma, lb = 0)
# estimate acceptance probability of crude Monte Carlo
print(attributes(est)$upbnd/est[1])
# reciprocal of acceptance probability
Z <- rtmvnorm(sigma = sigma, n = 1e3, lb = rep(0, n))
# sample exactly from auxiliary distribution
beta <- rtmvnorm(n = nrow(Z), sigma = C) + Z %*% X %*% C
# simulate beta given Z and plot boxplots of marginals
boxplot(beta, ylab = expression(beta))
# output the posterior means
colMeans(beta)

## End(Not run)

Multivariate truncated Student distribution

Description

Density, distribution function and random generation for the multivariate truncated Student distribution with location vector mu, scale matrix sigma, lower truncation limit lb, upper truncation limit ub and degrees of freedom df.

Arguments

Details

The truncation limits can include infinite values. The Monte Carlo (type = "mc") uses a sample of size B, while the qausi Monte Carlo (type = "qmc") uses a pointset of size ceiling(n/12) and estimates the relative error using 12 independent randomized QMC estimators.

pmvt computes an estimate and a deterministic upper bound of the distribution function of the multivariate normal distribution. Infinite values for vectors u and l are accepted. The Monte Carlo method uses sample size n: the larger n, the smaller the relative error of the estimator.

Value

dtmvt gives the density, ptmvt gives the distribution function, rtmvt generate random deviates.

Usage

dtmvt(x, mu, sigma, df, lb, ub, type = c("mc", "qmc"), log = FALSE, B = 1e4)
ptmvt(q, mu, sigma, df, lb, ub, type = c("mc", "qmc"), log = FALSE, B = 1e4)
rtmvt(n, mu, sigma, df, lb, ub)
pmvt(mu, sigma, df, lb = -Inf, ub = Inf, type = c("mc", "qmc"), B = 1e4)

Author(s)

Leo Belzile, R port from Matlab code by Z. I. Botev

References

Z. I. Botev and P. L'Ecuyer (2015), Efficient probability estimation and simulation of the truncated multivariate Student-t distribution, Proceedings of the 2015 Winter Simulation Conference, pp. 380-391

Examples

d <- 4; lb <- rep(0, d)
mu <- runif(d)
sigma <- matrix(0.5, d, d) + diag(0.5, d)
samp <- rtmvt(n = 10, mu = mu, sigma = sigma, df = 2, lb = lb)
loglik <- dtmvt(samp, mu = mu, sigma = sigma, df = 2, lb = lb, log = TRUE)
cdf <- ptmvt(samp, mu = mu, sigma = sigma, df = 2, lb = lb, log = TRUE, type = "q")

Truncated univariate normal distribution

Description

The function provides efficient state-of-the-art random number generation of a vector of truncated univariate distribution of the same length as the lower bound vector. The function is vectorized and the vector of means mu and of standard deviations sd are recycled.

If mu or sd are not specified they assume the default values of 0 and 1, respectively.

Arguments

Value

vector or matrix of random variates (rtnorm) or of quantiles (ptnorm), depending on the input

Examples

rtnorm(n = 10, mu = 2, lb = 1:10, ub = 2:11, method = "fast")
qtnorm(runif(10), mu = 2, lb = 1:10, ub = 2:11, sd = 1)

Fast truncated normal generator

Description

Efficient state-of-the-art generator of a vector of length(l)=length(u) from the standard multivariate normal distribution truncated over the region [l,u]. Infinite values for u and l are accepted.

Usage

trandn(l, u)

Arguments

Details

Suppose we wish to simulate a random variable Z drawn from N(\mu,\sigma^2) and conditional on l<Z<u using the inverse transform method. To achieve this, first compute X=norminvp(runif(1),(l-mu)/sig,(u-mu)/sig) and then set Z=mu+sig*X

Value

random variable drawn from the truncated normal distribution

Note

Use norminvp for the (slower) inverse transform method of simulating truncated normal variables.

Author(s)

Zdravko I. Botev

References

Z. I. Botev (2017), The Normal Law Under Linear Restrictions: Simulation and Estimation via Minimax Tilting, Journal of the Royal Statistical Society, Series B, 79 (1), pp. 1–24.

See Also

norminvp

Examples

trandn(l = 1,u = Inf)
trandn(l = rep(1, 10), u = rep(Inf, 10))

Truncated student generator for Bayesian regression simulation

Description

Simulates n random vectors X exactly distributed from the d-dimensional Student distribution with df=\nu degrees of freedom, mean zero and scale matrix sigma, conditional on l<X<u,

Usage

tregress(n, lb, ub, sigma, df)

Arguments

Value

list with components

• R: n vector of scale

• Z: a d by n matrix

so that \sqrt(\nu)Z/R follows a truncated Student distribution

Author(s)

Matlab code by Zdravko Botev, R port by Leo Belzile

References

Z. I. Botev and P. L'Ecuyer (2015), Efficient probability estimation and simulation of the truncated multivariate Student-t distribution, Proceedings of the 2015 Winter Simulation Conference, pp. 380-391,

Examples

d <- 5
tregress(lb =rep(-2, d), ub = rep(2, d), df = 3, n = 10,
  sigma = diag(0.5, d) + matrix(1, d, d))