
R Under development (unstable) (2026-03-12 r89611) -- "Unsuffered Consequences"
Copyright (C) 2026 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu

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")
> 
> pmvnorm <- function(lower = -Inf, upper = Inf, mean = rep(0, length(lower)), 
+                     corr = NULL, sigma = NULL, algorithm = GenzBretz(), ...) {
+ 
+     if (!inherits(algorithm, "GenzBretz"))
+         return(mvtnorm::pmvnorm(lower = lower, upper = upper, mean = mean, 
+                          sigma = sigma, corr = corr, algorithm = algorithm, ...))
+ 
+     args <- mvtnorm:::checkmvArgs(lower = lower, upper = upper, mean = mean, 
+                                   sigma = sigma, corr = corr)
+     if (args$uni)
+         return(mvtnorm::pmvnorm(lower = lower, upper = upper, mean = mean, 
+                          sigma = sigma, corr = corr, algorithm = algorithm, ...))
+ 
+     if (!is.null(args$corr)) args$sigma <- args$corr
+ 
+     Chol <- try(chol(args$sigma), silent = TRUE)
+     if (inherits(Chol, "try-error")) 
+        return(mvtnorm::pmvnorm(lower = lower, upper = upper, mean = mean, 
+                          sigma = sigma, corr = corr, algorithm = algorithm, ...))
+     Chol <- matrix(t(Chol)[lower.tri(Chol, diag = TRUE)], ncol = 1)
+     Chol <- ltMatrices(Chol, diag = TRUE, byrow = FALSE)
+ 
+     args$chol <- Chol
+     
+     M <- algorithm$maxpts
+     if (require("qrng", quietly = TRUE)) {
+         w <- t(ghalton(M, d = length(args$lower) - 1))
+     } else {
+         w <- NULL
+     }
+     args$w <- w
+     args$M <- M
+     args$seed <- 290875
+     args$logLik <- FALSE
+     args$corr <- args$sigma <- args$uni <- NULL
+     args$lower <- matrix(args$lower, ncol = 1)
+     args$upper <- matrix(args$upper, ncol = 1)
+ 
+     ret <- exp(do.call("lpmvnorm", args))
+     ret[ret < .Machine$double.eps] <- 0
+     ret
+ }
> 
> try(source("regtest-TVPACK.R", echo = TRUE))

> library("mvtnorm")

> chk <- function(...) isTRUE(all.equal(...))

> (cor1 <- toeplitz(c(1, 1/4, -1/8)))
       [,1] [,2]   [,3]
[1,]  1.000 0.25 -0.125
[2,]  0.250 1.00  0.250
[3,] -0.125 0.25  1.000

> (up1 <- c(1/4, 7/4, 5/8))
[1] 0.250 1.750 0.625

> d <- length(up1)

> pmvt.. <- function(df, algorithm) vapply(df, function(df) pmvt(upper = up1, 
+     corr = cor1, df = df, algorithm = algorithm), numeric(1))

> dfs <- 1:9

> pmvt_TV.7 <- replicate(7, pmvt..(dfs, TVPACK()))

> stopifnot(identical(unique(c(pmvt_TV.7 - pmvt_TV.7[, 
+     1])), 0))

> (pmvt.TV. <- pmvt_TV.7[, 1])
[1] 0.3554119 0.3817313 0.3923546 0.3980570 0.4016042 0.4040204 0.4057708
[8] 0.4070968 0.4081358

> (pmvt.TV <- pmvt..(dfs, TVPACK(1e-14)))
[1] 0.3554119 0.3817313 0.3923546 0.3980570 0.4016042 0.4040204 0.4057708
[8] 0.4070968 0.4081358

> chk(max(abs(pmvt.TV - pmvt.TV.)), 0)
[1] TRUE

> set.seed(47)

> pmvt_7 <- replicate(7, vapply(dfs, function(df) pmvt(df = df, 
+     upper = up1, corr = cor1), numeric(1)))

> relE <- 1 - pmvt_7/pmvt.TV

> rng.rE <- range(abs(relE))

> stopifnot(1e-06 < rng.rE[1], rng.rE[2] < 7e-04)

> stopifnot(chk(colMeans(abs(relE)), c(88, 64, 105, 
+     73, 52, 90, 87) * 1e-06, tol = 0.001))

> set.seed(29)

> corr <- cov2cor(crossprod(matrix(runif(9, -1, 1), 
+     3, 3)) + diag(3))

> df <- rpois(1, 3) + 1

> ctrl <- GenzBretz(maxpts = 2500000, abseps = 1e-06, 
+     releps = 0)

> upper <- rexp(3, 1)

> pmvt(upper = upper, corr = corr, df = df, algorithm = ctrl)
[1] 0.3920567
attr(,"error")
[1] 6.089669e-07
attr(,"msg")
[1] "Normal Completion"

> pmvt(upper = upper, corr = corr, df = df, algorithm = TVPACK())
[1] 0.3920566
attr(,"error")
[1] 1e-06
attr(,"msg")
[1] "Normal Completion"

> lower <- -rexp(3, 1)

> pmvt(lower = lower, upper = rep(Inf, 3), corr = corr, 
+     df = df, algorithm = ctrl)
[1] 0.4634843
attr(,"error")
[1] 6.46065e-07
attr(,"msg")
[1] "Normal Completion"

> pmvt(lower = lower, upper = rep(Inf, 3), corr = corr, 
+     df = df, algorithm = TVPACK())
[1] 0.4634844
attr(,"error")
[1] 1e-06
attr(,"msg")
[1] "Normal Completion"

> delt <- rexp(3, 1/10)

> upper <- delt + runif(3)

> ctrl <- GenzBretz(maxpts = 2500000, abseps = 1e-06, 
+     releps = 0)

> pmvt(upper = upper, corr = corr, df = df, algorithm = ctrl, 
+     delta = delt)
[1] 0.3235424
attr(,"error")
[1] 9.86407e-07
attr(,"msg")
[1] "Normal Completion"

> tools::assertError(pmvt(upper = upper, corr = corr, 
+     df = df, algorithm = TVPACK(), delta = delt))

> upper <- rexp(3, 1)

> pmvnorm(upper = upper, corr = corr, algorithm = ctrl)
[1] 0.7733946

> pmvnorm(upper = upper, corr = corr, algorithm = TVPACK())
[1] 0.7733949
attr(,"error")
[1] 1e-06
attr(,"msg")
[1] "Normal Completion"

> lower <- rexp(3, 5)

> pmvnorm(lower = lower, upper = rep(Inf, 3), corr = corr, 
+     algorithm = ctrl)
[1] 0.07507118

> pmvnorm(lower = lower, upper = rep(Inf, 3), corr = corr, 
+     algorithm = TVPACK())
[1] 0.07507119
attr(,"error")
[1] 1e-06
attr(,"msg")
[1] "Normal Completion"

> delt <- rexp(3, 1/10)

> upper <- delt + rexp(3, 1)

> pmvnorm(upper = upper, corr = corr, algorithm = ctrl, 
+     mean = delt)
[1] 0.3986757

> pmvnorm(upper = upper, corr = corr, algorithm = TVPACK(), 
+     mean = delt)
[1] 0.3986758
attr(,"error")
[1] 1e-06
attr(,"msg")
[1] "Normal Completion"

> corr <- cov2cor(crossprod(matrix(runif(4, -1, 1), 
+     2, 2)) + diag(2))

> upper <- rexp(2, 1)

> df <- rpois(1, runif(1, 0, 20))

> pmvt(upper = upper, corr = corr, df = df, algorithm = ctrl)
[1] 0.4565892
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"

> pmvt(upper = upper, corr = corr, df = df, algorithm = TVPACK())
[1] 0.4565892
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"

> pmvt(lower = -upper, upper = rep(Inf, 2), corr = corr, 
+     df = df, algorithm = ctrl)
[1] 0.4565892
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"

> pmvt(lower = -upper, upper = rep(Inf, 2), corr = corr, 
+     df = df, algorithm = TVPACK())
[1] 0.4565892
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"

> delt <- rexp(2, 1/5)

> upper <- delt + rexp(2, 1)

> pmvnorm(upper = upper, corr = corr, algorithm = ctrl, 
+     mean = delt)
[1] 0.320247

> pmvnorm(upper = upper, corr = corr, algorithm = TVPACK(), 
+     mean = delt)
[1] 0.320247
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"

> corr <- cov2cor(crossprod(matrix(runif(4, -1, 1), 
+     2, 2)) + diag(2))

> upper <- rexp(2, 1)

> pmvnorm(upper = upper, corr = corr, algorithm = Miwa(steps = 128))
[1] 0.4466655
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"

> pmvnorm(upper = upper, corr = corr, algorithm = TVPACK())
[1] 0.4466655
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"

> corr <- cov2cor(crossprod(matrix(runif(9, -1, 1), 
+     3, 3)) + diag(3))

> upper <- rexp(3, 1)

> ctrl <- Miwa(steps = 128)

> pmvnorm(upper = upper, corr = corr, algorithm = ctrl)
[1] 0.4499242
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"

> pmvnorm(upper = upper, corr = corr, algorithm = TVPACK())
[1] 0.4499242
attr(,"error")
[1] 1e-06
attr(,"msg")
[1] "Normal Completion"

> S <- toeplitz(c(1, 1/2, 1/4))

> set.seed(11)

> P0 <- pmvnorm(lower = c(-Inf, 0, 0), upper = Inf, 
+     corr = S)

> P1 <- pmvnorm(lower = c(-Inf, 0, 0), upper = Inf, 
+     corr = S, algorithm = TVPACK())

> P2 <- pmvnorm(lower = c(-Inf, 0, 0), upper = Inf, 
+     corr = S, algorithm = Miwa())

> P2a <- pmvnorm(lower = c(-Inf, 0, 0), upper = Inf, 
+     corr = S, algorithm = Miwa(512))

> P2. <- pmvnorm(lower = c(-Inf, 0, 0), upper = Inf, 
+     corr = S, algorithm = Miwa(2048))

> stopifnot(chk(1/3, c(P0), tol = 1e-14), chk(1/3, c(P1), 
+     tol = 1e-14), chk(1/3, c(P2), tol = 1e-09), chk(1/3, c(P2a), 
+     tol = 4e-12), chk .... [TRUNCATED] 
Error : chk(1/3, c(P0), tol = 1e-14) is not TRUE
> try(source("test-noisy-root.R", echo = TRUE))

> library("mvtnorm")

> set.seed(1)

> dim <- 10

> df <- 5

> D <- diag(dim) + crossprod(matrix(runif(25, -1, 1), 
+     dim, dim))

> corr <- cov2cor(D)

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

> 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.9500024
attr(,"error")
[1] 0.0001403885
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.9500143

> 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.9499675
attr(,"error")
[1] 0.0001052159
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.9500265

> 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.950041
attr(,"error")
[1] 0.000113758
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.9499789
> try(source("bugfix-tests.R", echo = TRUE))

> invisible(options(echo = TRUE))

> library("mvtnorm")

> set.seed(290875)

> chk <- function(...) isTRUE(all.equal(...))

> a <- 4.048

> shi <- -9

> slo <- -10

> mu <- -5

> sig <- matrix(c(1, 1, 1, 2), ncol = 2)

> pmvnorm(lower = c(-a, slo), upper = c(a, shi), mean = c(mu, 
+     2 * mu), sigma = sig)
[1] 0.04210637

> n <- 5

> lower <- -1

> upper <- 3

> df <- 4

> corr <- diag(5)

> corr[lower.tri(corr)] <- 0.5

> delta <- rep(0, 5)

> set.seed(290875)

> prob1 <- pmvt(lower = lower, upper = upper, delta = delta, 
+     df = df, corr = corr)

> set.seed(290875)

> prob2 <- pmvt(lower = lower, upper = upper, delta = delta, 
+     df = df, corr = corr)

> stopifnot(chk(prob1, prob2))

> a <- pmvnorm(lower = -Inf, upper = 2, mean = 0, sigma = matrix(1.5))

> attributes(a) <- NULL

> stopifnot(chk(a, pnorm(2, sd = sqrt(1.5))))

> a <- pmvnorm(lower = -Inf, upper = 2, mean = 0, sigma = matrix(0.5))

> attributes(a) <- NULL

> stopifnot(chk(a, pnorm(2, sd = sqrt(0.5))))

> a <- pmvnorm(lower = -Inf, upper = 2, mean = 0, sigma = 0.5)

> attributes(a) <- NULL

> stopifnot(chk(a, pnorm(2, sd = sqrt(0.5))))

> dmvnorm(x = c(0, 0), mean = c(1, 1), log = TRUE)
[1] -2.837877

> dmvnorm(x = c(0, 0), mean = c(25, 25), log = TRUE)
[1] -626.8379

> dmvnorm(x = c(0, 0), mean = c(30, 30), log = TRUE)
[1] -901.8379

> stopifnot(chk(dmvnorm(x = 0, mean = 30, log = TRUE), 
+     dnorm(0, 30, log = TRUE)))

> stopifnot(chk(f. <- dmvnorm(x = c(0, 0), mean = c(30, 
+     30), log = TRUE), dmvt(x = c(0, 0), delta = c(30, 30), log = TRUE, 
+     df = Inf)), c .... [TRUNCATED] 

> pnorm(2)^2
[1] 0.9550173

> pmvt(lower = c(-Inf, -Inf), upper = c(2, 2), delta = c(0, 
+     0), df = 25, corr = diag(2))
[1] 0.9446454
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"

> pmvt(lower = c(-Inf, -Inf), upper = c(2, 2), delta = c(0, 
+     0), df = 250, corr = diag(2))
[1] 0.9539846
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"

> pmvt(lower = c(-Inf, -Inf), upper = c(2, 2), delta = c(0, 
+     0), df = 1340, corr = diag(2))
[1] 0.9548248
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"

> pmvt(lower = c(-Inf, -Inf), upper = c(2, 2), delta = c(0, 
+     0), df = 2500, corr = diag(2))
[1] 0.9549141
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"

> pmvt(lower = c(-100, -100), upper = c(2, 2), delta = c(0, 
+     0), df = 2500, corr = diag(2))
[1] 0.9549141
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"

> pmvt(lower = c(-Inf, -Inf), upper = c(2, 2), delta = c(0, 
+     0), df = 0, corr = diag(2))
[1] 0.9550173
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"

> pmvt(lower = -Inf, upper = 2, delta = 0, df = 0, corr = 1)
    upper 
0.9772499 
attr(,"error")
[1] 0
attr(,"msg")
[1] "univariate: using pnorm"

> pnorm(2)
[1] 0.9772499

> pnorm(2)^2
[1] 0.9550173

> pmvnorm(lower = rep(-Inf, 2), upper = rep(2, 2), sigma = diag(2))
[1] 0.9550173

> pnorm(2)^90
[1] 0.1260393

> pmvnorm(lower = rep(-Inf, 90), upper = rep(2, 90), 
+     sigma = diag(90))
[1] 0.1260393

> pnorm(2)^199
[1] 0.01025932

> pmvnorm(lower = rep(-Inf, 199), upper = rep(2, 199), 
+     sigma = diag(199))
[1] 0.01025932

> cr = matrix(0.5, nr = 4, nc = 4)

> diag(cr) = 1

> cr
     [,1] [,2] [,3] [,4]
[1,]  1.0  0.5  0.5  0.5
[2,]  0.5  1.0  0.5  0.5
[3,]  0.5  0.5  1.0  0.5
[4,]  0.5  0.5  0.5  1.0

> a <- pmvt(low = -rep(1, 4), upp = rep(1, 4), df = 999, 
+     corr = cr)

> b <- pmvt(low = -rep(1, 4), upp = rep(1, 4), df = 4999, 
+     corr = cr)

> b
[1] 0.2892618
attr(,"error")
[1] 0.000110303
attr(,"msg")
[1] "Normal Completion"

> attributes(a) <- NULL

> attributes(b) <- NULL

> stopifnot(chk(round(a, 3), round(b, 3)))

> stopifnot(chk(c(pmvnorm(upper = c(-Inf, 1))), 0))

> stopifnot(chk(c(pmvnorm(lower = c(Inf, 1))), 0))

> stopifnot(chk(c(pmvnorm(lower = c(-2, 0), upper = c(-1, 
+     1), corr = matrix(rep(1, 4), 2, 2))), 0))

> stopifnot(chk(c(pmvnorm(-Inf, c(Inf, 0), 0, diag(2))), 
+     c(pmvnorm(-Inf, c(Inf, 0), 0))))

> stopifnot(chk(c(pmvnorm(lo = c(-Inf, -Inf), up = c(Inf, 
+     Inf), mean = c(0, 0))), 1))

> dm <- 250000

> iters <- 2

> corr <- 0.7

> dim <- 100

> abserr <- 3.5e-06

> cutoff <- -5.199338

> mn <- rep(0, dim)

> mat <- diag(dim)

> for (i in 1:dim) {
+     for (j in 1:(i - 1)) {
+         mat[i, j] = mat[j, i] = corr^(i - j)
+     }
+ }

> ll <- rep(cutoff, dim)

> mn <- rep(0, dim)

> p <- matrix(0, iters, 1)

> set.seed(290875)

> for (i in 1:iters) {
+     pp <- pmvnorm(lower = ll, sigma = mat, maxpts = dm, abseps = abserr)
+     p[i] <- 1 - pp
+ }

> stopifnot(abs(p[1] - p[2]) < 2 * abserr)

> ptmp <- p

> set.seed(290875)

> for (i in 1:iters) {
+     pp <- pmvnorm(lower = ll, sigma = mat, maxpts = dm, abseps = abserr)
+     p[i] <- 1 - pp
+ }

> stopifnot(chk(p, ptmp))

> pmvnormM <- function(...) pmvnorm(..., algorithm = Miwa())

> a <- 4.048

> shi <- -9

> slo <- -10

> mu <- -5

> sig <- matrix(c(1, 1, 1, 2), ncol = 2)

> pmvnormM(lower = c(-a, slo), upper = c(a, shi), mean = c(mu, 
+     2 * mu), sigma = sig)
[1] 0.04210555
attr(,"error")
[1] NA
attr(,"msg")
[1] "Normal Completion"

> n <- 5

> lower <- -1

> upper <- 3

> df <- 4

> corr <- diag(5)

> corr[lower.tri(corr)] <- 0.5

> delta <- rep(0, 5)

> set.seed(290875)

> prob1 <- pmvnormM(lower = lower, upper = upper, mean = delta, 
+     corr = corr)

> set.seed(290875)

> prob2 <- pmvnormM(lower = lower, upper = upper, mean = delta, 
+     corr = corr)

> stopifnot(chk(prob1, prob2))

> a <- pmvnormM(lower = -Inf, upper = 2, mean = 0, sigma = matrix(1.5))

> attributes(a) <- NULL

> stopifnot(chk(a, pnorm(2, sd = sqrt(1.5))))

> a <- pmvnormM(lower = -Inf, upper = 2, mean = 0, sigma = matrix(0.5))

> attributes(a) <- NULL

> stopifnot(chk(a, pnorm(2, sd = sqrt(0.5))))

> a <- pmvnormM(lower = -Inf, upper = 2, mean = 0, sigma = 0.5)

> attributes(a) <- NULL

> stopifnot(chk(a, pnorm(2, sd = sqrt(0.5))))

> stopifnot(chk(c(pmvnormM(upper = c(-Inf, 1))), 0))

> stopifnot(chk(c(pmvnormM(lower = c(Inf, 1))), 0))

> stopifnot(chk(pmvnormM(-Inf, c(Inf, 0), 0, diag(2)), 
+     pmvnormM(-Inf, c(Inf, 0), 0)))

> stopifnot(chk(c(pmvnormM(lo = c(-Inf, -Inf), up = c(Inf, 
+     Inf), mean = c(0, 0))), 1))

> dm <- 250000

> iters <- 2

> corr <- 0.7

> dim <- 10

> abserr <- 3.5e-06

> cutoff <- -5.199338

> mn <- rep(0, dim)

> mat <- diag(dim)

> for (i in 1:dim) {
+     for (j in 1:(i - 1)) {
+         mat[i, j] = mat[j, i] = corr^(i - j)
+     }
+ }

> ll <- rep(cutoff, dim)

> mn <- rep(0, dim)

> p <- matrix(0, iters, 1)

> set.seed(290875)

> for (i in 1:iters) {
+     pp <- pmvnormM(lower = ll, sigma = mat, maxpts = dm, abseps = abserr)
+     p[i] <- 1 - pp
+ }

> stopifnot(abs(p[1] - p[2]) < 2 * abserr)

> ptmp <- p

> set.seed(290875)

> for (i in 1:iters) {
+     pp <- pmvnormM(lower = ll, sigma = mat, maxpts = dm, abseps = abserr)
+     p[i] <- 1 - pp
+ }

> stopifnot(chk(p, ptmp))

> stopifnot(chk(c(pmvnorm(c(-Inf, -Inf, 0, 0))), 0.25))

> set.seed(290875)

> n <- 1e+05

> df <- rpois(1, 1/rexp(1, 1)) + 1

> dim <- rpois(1, runif(1, 0, 10)) + 2

> mn <- rnorm(dim, 0, 4)

> sigma <- matrix(runif(dim^2, -1, 1), nrow = dim, ncol = dim)

> sigma <- crossprod(sigma) + diag(dim)

> d <- runif(dim, 0.3, 20)

> sigma <- diag(d) %*% sigma %*% diag(d)

> corrMat <- cov2cor(sigma)

> sims1 <- rmvt(n, sigma = sigma, delta = mn, df = df, 
+     type = "shifted", pre0.9_9994 = TRUE)

> sims2 <- rmvt(n, sigma = sigma, delta = mn, df = df, 
+     type = "Kshirsagar", pre0.9_9994 = TRUE)

> lower <- mn - d * 2

> upper <- mn + d * 3

> comp <- function(x, lower, upper) {
+     all(x > lower) & all(x < upper)
+ }

> ind1 <- apply(sims1, 1, comp, lower = lower, upper = upper)

> mean(ind1)
[1] 0.247

> pmvt(lower, upper, sigma = sigma, delta = mn, df = df, 
+     type = "shifted")
[1] 0.2459638
attr(,"error")
[1] 8.791565e-05
attr(,"msg")
[1] "Normal Completion"

> ind2 <- apply(sims2, 1, comp, lower = lower, upper = upper)

> mean(ind2)
[1] 0.24729

> pmvt(lower, upper, sigma = sigma, delta = mn, df = df, 
+     type = "Kshirsagar")
[1] 0.2462984
attr(,"error")
[1] 6.430839e-05
attr(,"msg")
[1] "Normal Completion"

> sims1 <- rmvt(n, sigma = corrMat, delta = mn, df = df, 
+     type = "shifted", pre0.9_9994 = TRUE)

> sims2 <- rmvt(n, sigma = corrMat, delta = mn, df = df, 
+     type = "Kshirsagar", pre0.9_9994 = TRUE)

> lower <- mn - d * 0.5

> upper <- mn + d

> comp <- function(x, lower, upper) {
+     all(x > lower) & all(x < upper)
+ }

> ind1 <- apply(sims1, 1, comp, lower = lower, upper = upper)

> mean(ind1)
[1] 0.99669

> pmvt(lower, upper, corr = corrMat, delta = mn, df = df, 
+     type = "shifted")
[1] 0.996872
attr(,"error")
[1] 2.461945e-05
attr(,"msg")
[1] "Normal Completion"

> ind2 <- apply(sims2, 1, comp, lower = lower, upper = upper)

> mean(ind2)
[1] 0.98827

> pmvt(lower, upper, corr = corrMat, delta = mn, df = df, 
+     type = "Kshirsagar")
[1] 0.9882905
attr(,"error")
[1] 6.344011e-05
attr(,"msg")
[1] "Normal Completion"

> m <- 10

> rho <- 0.1

> k <- 2

> alpha <- 0.05

> cc_z <- numeric(m)

> var <- matrix(c(1, rho, rho, 1), nrow = 2, ncol = 2, 
+     byrow = T)

> i <- 1

> q1 <- qmvnorm((k * (k - 1))/(m * (m - 1)) * alpha, 
+     tail = "upper", sigma = var, ptol = 1e-05)$quantile

> q2 <- qmvnorm((k * (k - 1))/(m * (m - 1)) * alpha, 
+     tail = "upper", sigma = var, interval = c(0, 5), ptol = 1e-05)$quantile

> stopifnot(chk(round(q1, 4), round(q2, 4)))

> qmvnorm(0.95, sigma = tcrossprod(c(0.009, 0.75, 0.25)))$quantile
[1] 1.23364

> stopifnot(is.finite(qmvt(0.95, df = 0, corr = matrix(1))$quantile))

> corr <- matrix(1, ncol = 2, nrow = 2)

> p <- c(pmvnorm(lower = c(-Inf, -Inf), upper = c(1.96, 
+     1.96), mean = c(1.72, 1.72), corr = corr), pmvt(lower = c(-Inf, 
+     -Inf), upper = c .... [TRUNCATED] 

> stopifnot(all(abs(p - mean(p)) < 1/100))

> m <- 3

> S <- diag(m)

> S[2, 1] <- S[1, 2] <- 1/4

> S[3, 1] <- S[1, 3] <- 1/5

> S[3, 2] <- S[2, 3] <- 1/3

> p <- pmvnorm(lower = c(-Inf, 0, 0), upper = c(0, Inf, 
+     Inf), mean = c(0, 0, 0), sigma = S, algorithm = Miwa())

> stopifnot(!is.na(p))

> set.seed(29)

> d1 <- function(x, mean, sigma) {
+     distval <- mahalanobis(x, center = mean, cov = sigma)
+     logdet <- sum(log(eigen(sigma, symmetric = TRUE,  .... [TRUNCATED] 

> d2 <- function(...) dmvnorm(..., log = TRUE)

> for (i in 1:100) {
+     p <- sample(2:10, 1)
+     Sigma <- tcrossprod(matrix(runif(p^2) * 2, ncol = p))
+     x <- matrix(rnorm(p), nr = 1)
+      .... [TRUNCATED] 

> L <- diag(10 * (1:4))

> L[lower.tri(L)] <- 1:6

> L[3, 3] <- 0

> L
     [,1] [,2] [,3] [,4]
[1,]   10    0    0    0
[2,]    1   20    0    0
[3,]    2    4    0    0
[4,]    3    5    6   40

> (Sig <- tcrossprod(L))
     [,1] [,2] [,3] [,4]
[1,]  100   10   20   30
[2,]   10  401   82  103
[3,]   20   82   20   26
[4,]   30  103   26 1670

> set.seed(123)

> fx <- dmvnorm(rbind(0, 1:4, matrix(rnorm(4 * 10), 
+     ncol = 4)), sigma = Sig)

> stopifnot(chk(fx, c(Inf, rep(0, 1 + 10))))

> ret <- structure(list(N = 10, NU = 25, LOWER = c(-0.430060315238938, 
+     -0.430060315238938, -0.430060315238938, -0.430060315238938, 
+     -0.43 .... [TRUNCATED] 

> RS <- c(403, 480, 641015092, 1848202935, -2124158291, 
+     -2116162620, 1818211306, -796165035, -1592745489, -483415562, 
+     -77025504, -170853 .... [TRUNCATED] 

> f <- function() {
+     error <- 0
+     value <- 0
+     inform <- 0
+     ret <- .C(C_mvtdst, N = as.integer(ret$N), NU = as.integer(ret$NU), 
+   .... [TRUNCATED] 

> .Random.seed <- RS

> p <- 0.95

> stopifnot(chk(round(qmvnorm(p, sigma = diag(3), tail = "upper")$quantile, 
+     2), round(qnorm(p^(1/3), lower = FALSE), 2)))

> stopifnot(chk(round(qmvnorm(p, sigma = diag(3), tail = "lower")$quantile, 
+     2), round(qnorm(p^(1/3), lower = TRUE), 2)))

> set.seed(29)

> p <- 0.95

> d <- 4

> qmvnorm(p, sigma = diag(d), tail = "lower")$quantile
[1] 2.23395

> qmvnorm(p, sigma = diag(d), tail = "upper")$quantile
[1] -2.23395

> qmvnorm(p, sigma = diag(d), tail = "both")$quantile
[1] 2.490844

> p <- 1 - 0.95

> d <- 4

> qmvnorm(p, sigma = diag(d), tail = "lower")$quantile
[1] -0.06784195

> qmvnorm(p, sigma = diag(d), tail = "upper")$quantile
[1] 0.06784195

> qmvnorm(p = 0.5, tail = "lower", mean = c(6.75044368, 
+     0.04996326), sigma = rbind(c(0.1026055, 0.02096418), c(0.02096418, 
+     0.16049956))) .... [TRUNCATED] 
[1] 6.750444

> stint <- c(6.75044332319072, 6.75044368)

> qmvnorm(p = 0.5, tail = "lower", mean = c(6.75044368, 
+     0.04996326), sigma = rbind(c(0.1026055, 0.02096418), c(0.02096418, 
+     0.16049956)), .... [TRUNCATED] 
[1] 6.750444

> R2 = matrix(c(0.7071068, 0.6924398, 0.7054602, 0.7054602, 
+     0.6292745, 0.6924398, 0.7071068, 0.6909812, 0.6909712, 0.612867, 
+     0.7054602,  .... [TRUNCATED] 

> call <- try(qmvnorm(p = 1 - 0.0001726701, mean = c(-0.8752332, 
+     -0.9487915, -0.9719237, -0.5855204, -0.9046457), sigma = R2, 
+     tail = "lo ..." ... [TRUNCATED] 

> inherits(call, "try-error")
[1] TRUE

> grepl("Covariance matrix not positive semidefinite", 
+     geterrmessage())
[1] TRUE

> call <- try(qmvt(p = 1 - 0.0001726701, mean = c(-0.8752332, 
+     -0.9487915, -0.9719237, -0.5855204, -0.9046457), sigma = R2, 
+     tail = "lower ..." ... [TRUNCATED] 

> inherits(call, "try-error")
[1] TRUE

> grepl("Covariance matrix not positive semidefinite", 
+     geterrmessage())
[1] TRUE

> all.equal(qnorm(p = 0.2397501, mean = 1, sd = sqrt(2)), 
+     qmvnorm(p = 0.2397501, mean = 1, sigma = 2)$quantile)
[1] TRUE
> 
> proc.time()
   user  system elapsed 
 17.255   0.075  17.334 
