| Title: | Sampling from Truncated Multivariate Normal and t Distributions |
|---|---|
| Description: | Efficient sampling of truncated multivariate (scale) mixtures of normals under linear inequality constraints is nontrivial due to the analytically intractable normalizing constant. Meanwhile, traditional methods may subject to numerical issues, especially when the dimension is high and dependence is strong. Algorithms proposed by Li and Ghosh (2015) <doi: 10.1080/15598608.2014.996690> are adopted for overcoming difficulties in simulating truncated distributions. Efficient rejection sampling for simulating truncated univariate normal distribution is included in the package, which shows superiority in terms of acceptance rate and numerical stability compared to existing methods and R packages. An efficient function for sampling from truncated multivariate normal distribution subject to convex polytope restriction regions based on Gibbs sampler for conditional truncated univariate distribution is provided. By extending the sampling method, a function for sampling truncated multivariate Student's t distribution is also developed. Moreover, the proposed method and computation remain valid for high dimensional and strong dependence scenarios. Empirical results in Li and Ghosh (2015) <doi: 10.1080/15598608.2014.996690> illustrated the superior performance in terms of various criteria (e.g. mixing and integrated auto-correlation time). |
| Authors: | Ting Fung (Ralph) Ma [cre, aut], Sujit K. Ghosh [aut], Yifang Li [aut] |
| Maintainer: | Ting Fung (Ralph) Ma <[email protected]> |
| License: | GPL-2 |
| Version: | 1.1.1 |
| Built: | 2024-12-01 08:39:33 UTC |
| Source: | CRAN |
dtuvn calculates the probability density function (pdf) of truncated univariate normal distribution.
dtuvn(x, mean, sd, lower, upper)dtuvn(x, mean, sd, lower, upper)
x |
value at which density is desired. |
mean |
mean of the underlying univariate normal distribution. |
sd |
standard deviation of the underlying univariate normal distribution. |
lower |
lower bound for truncation. |
upper |
upper bound for truncation. |
dtuvn returns the density (with same dimension and type as x) of truncated univariate normal distribution.
dtuvn(x= -3:3, mean=0, sd=1 ,lower= -2, upper=2)dtuvn(x= -3:3, mean=0, sd=1 ,lower= -2, upper=2)
exp_acc_opt calculates the acceptance rate of translated-exponential rejection sampling for the truncation interval (a,b).
exp_acc_opt(a, b)exp_acc_opt(a, b)
a |
lower bound for truncation. |
b |
upper bound for truncation. |
set.seed(1203) exp_acc_opt(1,2)set.seed(1203) exp_acc_opt(1,2)
exp_rej is used for translated-exponential rejection sampling.
exp_rej(a, b = Inf, lam = "default")exp_rej(a, b = Inf, lam = "default")
a |
lower bound |
b |
upper bound |
lam |
lambda for translated-exponential only |
exp_rej returns a list
x: sampled value; and
acc: total number of draw used.
set.seed(1) exp_rej(a=1, b=Inf)set.seed(1) exp_rej(a=1, b=Inf)
halfnorm_acc calculates the acceptance rate of half-normal rejection sampling for the truncation interval (a,b).
halfnorm_acc(a, b)halfnorm_acc(a, b)
a |
lower bound for truncation. |
b |
upper bound for truncation. |
set.seed(1203) halfnorm_acc(1,2)set.seed(1203) halfnorm_acc(1,2)
halfnorm_rej is used for half-normal rejection sampling.
halfnorm_rej(a, b)halfnorm_rej(a, b)
a |
lower bound |
b |
upper bound |
halfnorm_rej returns a list
x: sampled value; and
acc: total number of draw used.
set.seed(1) halfnorm_rej(a=1, b=Inf)set.seed(1) halfnorm_rej(a=1, b=Inf)
imp contains a general function for rejection sampling of standardized truncated univariate normal distribution in (a,b).
imp(a, b)imp(a, b)
a |
lower bound for truncation. |
b |
upper bound for truncation (must be |
imp returns a list
x: sampled value; and
acc: total number of draw used.
imp(1,Inf) # Case 1: [a,infty) imp(-1,1) # Case 2: 0 in [a,b], a<0<b imp(1,2) # Case 3: [a,b], a>=0imp(1,Inf) # Case 1: [a,infty) imp(-1,1) # Case 2: 0 in [a,b], a<0<b imp(1,2) # Case 3: [a,b], a>=0
imp_acc calculates the acceptance rate of truncated univariate standardized normal distribution rejection sampling for the truncation interval (a,b).
imp_acc(a, b)imp_acc(a, b)
a |
lower bound for truncation. |
b |
upper bound for truncation. |
imp_acc(1,Inf) # Case 1: [a,infty) imp_acc(-1,1) # Case 2: 0 in [a,b], a<0<b imp_acc(1,2) # Case 3: [a,b], a>=0imp_acc(1,Inf) # Case 1: [a,infty) imp_acc(-1,1) # Case 2: 0 in [a,b], a<0<b imp_acc(1,2) # Case 3: [a,b], a>=0
norm_acc calculates the acceptance rate of normal rejection sampling for the truncation interval (a,b).
norm_acc(a, b)norm_acc(a, b)
a |
lower bound for truncation. |
b |
upper bound for truncation. |
set.seed(1203) norm_acc(1,2)set.seed(1203) norm_acc(1,2)
norm_rej is used for normal rejection sampling.
norm_rej(a, b = Inf)norm_rej(a, b = Inf)
a |
lower bound |
b |
upper bound |
norm_rej returns a list
x: sampled value; and
acc: total number of draw used.
set.seed(1) norm_rej(a=1, b=Inf)set.seed(1) norm_rej(a=1, b=Inf)
ptuvn calculates the cumulative distribution function (cdf) of truncated univariate normal distribution.
ptuvn(x, mean, sd, lower, upper)ptuvn(x, mean, sd, lower, upper)
x |
value at which cdf is desired. |
mean |
mean of the underlying univariate normal distribution. |
sd |
standard deviation of the underlying univariate normal distribution. |
lower |
lower bound for truncation. |
upper |
upper bound for truncation. |
ptuvn returns the cumulative distribution function (with same dimension and type as x) of truncated univariate normal distribution.
ptuvn(x= -3:3, mean=0, sd=1 ,lower= -2, upper=2)ptuvn(x= -3:3, mean=0, sd=1 ,lower= -2, upper=2)
rtmvn simulates truncated multivariate (p-dimensional) normal distribution subject to linear inequality constraints. The constraints should be written as a matrix (D) with lower and upper as the lower and upper bounds for those constraints respectively. Note that D can be non-full rank, which generalize many traditional methods.
rtmvn( n, Mean, Sigma, D = diag(1, length(Mean)), lower, upper, int = NULL, burn = 10, thin = 1 )rtmvn( n, Mean, Sigma, D = diag(1, length(Mean)), lower, upper, int = NULL, burn = 10, thin = 1 )
n |
number of random samples desired (sample size). |
Mean |
mean vector of the underlying multivariate normal distribution. |
Sigma |
positive definite covariance matrix of the underlying multivariate normal distribution. |
D |
matrix or vector of coefficients of linear inequality constraints. |
lower |
vector of lower bounds for truncation. |
upper |
vector of upper bounds for truncation. |
int |
initial value vector for Gibbs sampler (satisfying truncation), if |
burn |
burn-in iterations discarded (default as |
thin |
thinning lag (default as |
rtmvn returns a (n*p) matrix (or vector when n=1) containing random numbers which approximately follows truncated multivariate normal distribution.
# Example for full rank with strong dependence d <- 3 rho <- 0.9 Sigma <- matrix(0, nrow=d, ncol=d) Sigma <- rho^abs(row(Sigma) - col(Sigma)) D1 <- diag(1,d) # Full rank set.seed(1203) ans.1 <- rtmvn(n=1000, Mean=1:d, Sigma, D=D1, lower=rep(-1,d), upper=rep(1,d), int=rep(0,d), burn=50) apply(ans.1, 2, summary) # Example for non-full rank d <- 3 rho <- 0.5 Sigma <- matrix(0, nrow=d, ncol=d) Sigma <- rho^abs(row(Sigma) - col(Sigma)) D2 <- matrix(c(1,1,1,0,1,0,1,0,1),ncol=d) qr(D2)$rank # 2 set.seed(1228) ans.2 <- rtmvn(n=100, Mean=1:d, Sigma, D=D2, lower=rep(-1,d), upper=rep(1,d), burn=10) apply(ans.2, 2, summary)# Example for full rank with strong dependence d <- 3 rho <- 0.9 Sigma <- matrix(0, nrow=d, ncol=d) Sigma <- rho^abs(row(Sigma) - col(Sigma)) D1 <- diag(1,d) # Full rank set.seed(1203) ans.1 <- rtmvn(n=1000, Mean=1:d, Sigma, D=D1, lower=rep(-1,d), upper=rep(1,d), int=rep(0,d), burn=50) apply(ans.1, 2, summary) # Example for non-full rank d <- 3 rho <- 0.5 Sigma <- matrix(0, nrow=d, ncol=d) Sigma <- rho^abs(row(Sigma) - col(Sigma)) D2 <- matrix(c(1,1,1,0,1,0,1,0,1),ncol=d) qr(D2)$rank # 2 set.seed(1228) ans.2 <- rtmvn(n=100, Mean=1:d, Sigma, D=D2, lower=rep(-1,d), upper=rep(1,d), burn=10) apply(ans.2, 2, summary)
rtmvt simulates truncated multivariate (p-dimensional) Student's t distribution subject to linear inequality constraints. The constraints should be written as a matrix (D) with lower and upper as the lower and upper bounds for those constraints respectively. Note that D can be non-full rank, which generalizes many traditional methods.
rtmvt(n, Mean, Sigma, nu, D, lower, upper, int = NULL, burn = 10, thin = 1)rtmvt(n, Mean, Sigma, nu, D, lower, upper, int = NULL, burn = 10, thin = 1)
n |
number of random samples desired (sample size). |
Mean |
location vector of the multivariate Student's t distribution. |
Sigma |
positive definite dispersion matrix of the multivariate t distribution. |
nu |
degrees of freedom for Student-t distribution. |
D |
matrix or vector of coefficients of linear inequality constraints. |
lower |
lower bound vector for truncation. |
upper |
upper bound vector for truncation. |
int |
initial value vector for Gibbs sampler (satisfying truncation), if |
burn |
burn-in iterations discarded (default as |
thin |
thinning lag (default as |
rtmvt returns a (n*p) matrix (or vector when n=1) containing random numbers which follows truncated multivariate Student-t distribution.
# Example for full rank d <- 3 rho <- 0.5 nu <- 10 Sigma <- matrix(0, nrow=d, ncol=d) Sigma <- rho^abs(row(Sigma) - col(Sigma)) D1 <- diag(1,d) # Full rank set.seed(1203) ans.t <- rtmvt(n=1000, Mean=1:d, Sigma, nu=nu, D=D1, lower=rep(-1,d), upper=rep(1,d), burn=50, thin=0) apply(ans.t, 2, summary)# Example for full rank d <- 3 rho <- 0.5 nu <- 10 Sigma <- matrix(0, nrow=d, ncol=d) Sigma <- rho^abs(row(Sigma) - col(Sigma)) D1 <- diag(1,d) # Full rank set.seed(1203) ans.t <- rtmvt(n=1000, Mean=1:d, Sigma, nu=nu, D=D1, lower=rep(-1,d), upper=rep(1,d), burn=50, thin=0) apply(ans.t, 2, summary)
rtuvn simulates truncated univariate normal distribution within the interval.
rtuvn(n = 1, mean = 0, sd = 1, lower, upper)rtuvn(n = 1, mean = 0, sd = 1, lower, upper)
n |
number of random samples desired (sample size). |
mean |
mean of the underlying univariate normal distribution. |
sd |
standard deviation of the underlying univariate normal distribution. |
lower |
lower bound for truncation. |
upper |
upper bound for truncation. |
rtuvn returns a vector of random number follows truncated univariate normal distribution.
set.seed(1203) ans <- rtuvn(n=1000, mean=1, sd=2, lower=-2, upper=3) summary(ans) # Check if the sample matches with CDF by KS test ks.test(ans,"ptuvn",1,2,-2,3)set.seed(1203) ans <- rtuvn(n=1000, mean=1, sd=2, lower=-2, upper=3) summary(ans) # Check if the sample matches with CDF by KS test ks.test(ans,"ptuvn",1,2,-2,3)
unif_acc calculates the acceptance rate of uniform rejection sampling for the truncation interval (a,b).
unif_acc(a, b)unif_acc(a, b)
a |
lower bound for truncation. |
b |
upper bound for truncation. |
set.seed(1203) unif_acc(1,2)set.seed(1203) unif_acc(1,2)
unif_rej is used for uniform rejection sampling.
unif_rej(a, b)unif_rej(a, b)
a |
lower bound |
b |
upper bound |
unif_rej returns a list
x: sampled value; and
acc: total number of draw used.
set.seed(1) unif_rej(a=1, b=2)set.seed(1) unif_rej(a=1, b=2)