R Under development (unstable) (2023-11-26 r85638 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 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") > > 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 Warning messages: 1: In probval.Miwa(algorithm, n, df, lower, upper, infin, corr, delta) : Approximating +/-Inf by +/-1000 2: In probval.Miwa(algorithm, n, df, lower, upper, infin, corr, delta) : Approximating +/-Inf by +/-1000 > > proc.time() user system elapsed 42.43 0.68 43.10