R Under development (unstable) (2024-08-16 r87026 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > > library("mvtnorm") > > set.seed(1) > dim <- 10 > df <- 5 > > D <- diag(dim)+crossprod(matrix(runif(25,-1,1),dim,dim)) > corr <- cov2cor(D) > > ## one-sided, lower tail > qu <- qmvt(0.95, df=df, corr=corr)$quantile > pmvt(lower=rep(-Inf, dim), upper=rep(qu,dim), corr=corr, df=df) [1] 0.9499886 attr(,"error") [1] 8.394993e-05 attr(,"msg") [1] "Normal Completion" > > qu <- qmvnorm(0.95, corr=corr)$quantile > pmvnorm(lower=rep(-Inf, dim), upper=rep(qu,dim), corr=corr) [1] 0.9499785 attr(,"error") [1] 6.265055e-05 attr(,"msg") [1] "Normal Completion" > > ## two-sided > qu <- qmvt(0.95, df=df, corr=corr, tail="both.tails")$quantile > pmvt(lower=rep(-qu, dim), upper=rep(qu,dim), corr=corr, df=df) [1] 0.949925 attr(,"error") [1] 9.635837e-05 attr(,"msg") [1] "Normal Completion" > > qu <- qmvnorm(0.95, corr=corr, tail="both.tails")$quantile > pmvnorm(lower=rep(-qu, dim), upper=rep(qu,dim), corr=corr) [1] 0.949927 attr(,"error") [1] 0.0001545113 attr(,"msg") [1] "Normal Completion" > > ## one-sided, upper tail > qu <- qmvt(0.95, df=df, corr=corr, tail="upper.tail")$quantile > pmvt(lower=rep(qu, dim), upper=rep(Inf,dim), corr=corr, df=df) [1] 0.9499528 attr(,"error") [1] 9.714629e-05 attr(,"msg") [1] "Normal Completion" > > qu <- qmvnorm(0.95, corr=corr, tail="upper.tail")$quantile > pmvnorm(lower=rep(qu, dim), upper=rep(Inf,dim), corr=corr) [1] 0.9500897 attr(,"error") [1] 7.617415e-05 attr(,"msg") [1] "Normal Completion" > > ## cross-check interval works > ## one-sided, lower tail > qu <- qmvt(0.95, df=df , corr=corr, interval=c(0,10))$quantile > pmvt(lower=rep(-Inf, dim), upper=rep(qu,dim), corr=corr, df=df) [1] 0.9500624 attr(,"error") [1] 0.0001741215 attr(,"msg") [1] "Normal Completion" > > qu <- qmvnorm(0.95, corr=corr, interval=c(0,10))$quantile > pmvnorm(lower=rep(-Inf, dim), upper=rep(qu,dim), corr=corr) [1] 0.9500022 attr(,"error") [1] 9.231686e-05 attr(,"msg") [1] "Normal Completion" > > proc.time() user system elapsed 8.85 0.09 8.93