R version 4.2.2 (2022-10-31) -- "Innocent and Trusting"
Copyright (C) 2022 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

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.

> invisible(options(echo = TRUE))
> library("mvtnorm")
> set.seed(290875)
> 
> chk <- function(...) isTRUE(all.equal(...))
> 
> # correlation matrices for unequal variances were wrong
> # from Pamela Ohman-Strickland <ohmanpa@UMDNJ.EDU>
> 
> 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.04210555
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"
> 
> # check if set.seed works (starting from 0.5-7)
> 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))
> 
> # confusion for univariate probabilities when sigma is a matrix
> # by Jerome Asselin <jerome@hivnet.ubc.ca>
> 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(.5))
> attributes(a) <- NULL
> stopifnot(chk(a, pnorm(2, sd=sqrt(.5))))
> a <- pmvnorm(lower=-Inf,upper=2,mean=0,sigma=.5)
> attributes(a) <- NULL
> stopifnot(chk(a, pnorm(2, sd=sqrt(.5))))
> 
> # log argument added by Jerome Asselin <jerome@hivnet.ubc.ca>
> 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(dmvnorm(x=c(0,0), mean =c(30,30),log=TRUE) -> f.,
+               dmvt   (x=c(0,0), delta=c(30,30),log=TRUE, df=Inf))
+     ,
+     chk(f., dmvt(x=c(0,0), delta=c(30,30),log=TRUE, df=10000),
+               tolerance = 0.09)
+ )
> 
> # large df
> 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"
> 
> # df = 0
> 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
> 
> # larger dimensions
> pnorm(2)^2
[1] 0.9550173
> pmvnorm(lower=rep(-Inf, 2), upper=rep(2,2), sigma = diag(2))
[1] 0.9550173
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"
> pnorm(2)^90
[1] 0.1260393
> pmvnorm(lower=rep(-Inf, 90), upper=rep(2,90), sigma = diag(90))
[1] 0.1260393
attr(,"error")
[1] 0
attr(,"msg")
[1] "Normal Completion"
> pnorm(2)^199
[1] 0.01025932
> pmvnorm(lower=rep(-Inf, 199), upper=rep(2,199), sigma = diag(199))
[1] 0.01025932
attr(,"error")
[1] 0
attr(,"msg")
[1] "Normal Completion"
> 
> # larger dimensions, again. Spotted by Chihiro Kuroki <kuroki@oak.dti.ne.jp>
> # Alan's fix to MVCHNC solves this problem
> 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.2893192
attr(,"error")
[1] 8.971803e-05
attr(,"msg")
[1] "Normal Completion"
> attributes(a) <- NULL
> attributes(b) <- NULL
> stopifnot(chk(round(a, 3), round(b, 3)))
> 
> # cases where the support is the empty set tried to compute something.
> # spotted by Peter Thomson <peter@statsresearch.co.nz>
> 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))
> 
> # bugged Fritz (long time ago)
> stopifnot(chk(c(pmvnorm(-Inf, c(Inf, 0), 0, diag(2))),
+ 		    c(pmvnorm(-Inf, c(Inf, 0), 0))))
> 
> # this is a bug in `mvtdst' nobody was able to fix yet :-(
> stopifnot(chk(c(pmvnorm(lo=c(-Inf,-Inf), up=c(Inf,Inf), mean=c(0,0))), 1))
> 
> ### check for correct random seed initialization
> ### problem reported by Karen Conneely <conneely@umich.edu>
> dm <- 250000
> iters <- 2
> corr <- .7
> dim <- 100
> abserr <- .0000035
> 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))
> 
> ### same for algoritm = Miwa
> 
> 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"
> 
> # check if set.seed works (starting from 0.5-7)
> 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))
> 
> # confusion for univariate probabilities when sigma is a matrix
> # by Jerome Asselin <jerome@hivnet.ubc.ca>
> 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(.5))
> attributes(a) <- NULL
> stopifnot(chk(a, pnorm(2, sd=sqrt(.5))))
> a <- pmvnormM(lower=-Inf,upper=2,mean=0,sigma=.5)
> attributes(a) <- NULL
> stopifnot(chk(a, pnorm(2, sd=sqrt(.5))))
> 
> 
> # cases where the support is the empty set tried to compute something.
> # spotted by Peter Thomson <peter@statsresearch.co.nz>
> stopifnot(chk(c(pmvnormM(upper=c(-Inf,1))), 0))
> stopifnot(chk(c(pmvnormM(lower=c(Inf,1))), 0))
> 
> # bugged Fritz (long time ago)
> stopifnot(chk(pmvnormM(-Inf, c(Inf, 0), 0, diag(2)),
+ 		    pmvnormM(-Inf, c(Inf, 0), 0)))
> 
> # this is a bug in `mvtdst' nobody was able to fix yet :-(
> stopifnot(chk(c(pmvnormM(lo=c(-Inf,-Inf), up=c(Inf,Inf), mean=c(0,0))), 1))
> 
> ### check for correct random seed initialization
> ### problem reported by Karen Conneely <conneely@umich.edu>
> dm <- 250000
> iters <- 2
> corr <- .7
> dim <- 10
> abserr <- .0000035
> 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))
> 
> ### was == 1; spotted by Alex Lenkoski <lenkoski@stat.washington.edu>
> stopifnot(chk(c(pmvnorm(c(-Inf, -Inf, 0, 0))), 0.25))
> 
> #############################
> ## testing rmvt und pmvt
> #############################
> set.seed(290875)
> n <- 100000
> df <- rpois(1,1/rexp(1,1))+1
> dim <- rpois(1,runif(1,0,10))+2
> mn <- rnorm(dim,0,4) ##rep(0,dim)
> 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)
> 
> ## sigma handed over
> 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) #Monte Carlo Integration
[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"
> 
> ## corrMat handed over
> 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) #Monte Carlo Integration
[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"
> 
> ### approx_interval for tail = "upper" went wild
> ### spotted by Ravi Varadhan <rvaradhan@jhmi.edu>
> 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=0.00001)$quantile
> q2 <- qmvnorm((k*(k-1))/(m*(m-1))*alpha, tail="upper", sigma=var,
+          interval = c(0, 5), ptol=0.00001)$quantile
> stopifnot(chk(round(q1, 4), round(q2, 4)))
> 
> ### grrr, still problems in approx_interval
> qmvnorm(.95, sigma = tcrossprod(c(0.009, 0.75, 0.25)))$quantile
[1] 1.23364
> 
> ### qmvt(..., df = 0, ...) didn't work
> ### spotted by Ulrich Halekoh <Ulrich.Halekoh@agrsci.dk>
> stopifnot(is.finite(qmvt(.95, df = 0, corr = matrix(1))$quantile))
> 
> ### spotted by <Tobias.Mielke@aptivsolutions.com> and fixed
> ### in mvtdst.f by Alan 2013-05-29
> 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(1.96,1.96),delta=c(1.72,1.72),df=0,corr=corr),
+        pmvt(lower=c(-Inf,-Inf),upper=c(1.96,1.96) - c(1.72,1.72),df=0,corr=corr),
+        pmvt(lower=c(-Inf,-Inf),upper=c(1.96,1.96) - c(1.72,1.72),df=100,corr=corr),
+        pmvt(lower=c(-Inf,-Inf),upper=c(1.96,1.96), delta=c(1.72,1.72),df=100,corr=corr))
> stopifnot(all(abs(p - mean(p)) < 1 / 100))
> 
> ### spotted and fixed by Xuefei Mi
> 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
> # NaN was given.
> 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))
> 
> ### introduced with dmvnorm in 0.9-9999
> set.seed(29)
> ### dmvnorm up to 0.9-9997
> d1 <- function(x, mean, sigma) {
+     distval <- mahalanobis(x, center = mean, cov = sigma)
+     logdet <- sum(log(eigen(sigma, symmetric=TRUE,
+                                    only.values=TRUE)$values))
+     -(ncol(x)*log(2*pi) + logdet + distval)/2
+ }
> ### current version
> 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)
+   m <- runif(p)
+   ld1 <- d1(x=x, mean=m, sigma=Sigma)
+   ld2 <- d2(x=x, mean=m, sigma=Sigma)
+ 
+   stopifnot(chk(ld1, ld2, tol = .Machine$double.eps^(1/3)))
+ }
> 
> ### --- Singular Sigma --- Now treated the same as  dnorm(*, sd=0):  "Inf or 0"
> L <- diag(10*(1:4))
> L[lower.tri(L)] <- 1:6
> L[3,3] <- 0 # to make it singular
> 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))))
> ## gave NaN for a long time, then error, then NaN, now we have converged ;-)
> 
> ### NaN spotted by David Charles Airey <airey_david_charles@lilly.com>
> ### data contains all input parameters and a special seed
> ret <-
+ structure(list(N = 10L, NU = 25L, LOWER = c(-0.430060315238938, 
+ -0.430060315238938, -0.430060315238938, -0.430060315238938, -0.430060315238938, 
+ -0.430060315238938, -0.430060315238938, -0.430060315238938, -0.430060315238938, 
+ -0.430060315238938), UPPER = c(0.430060315238938, 0.430060315238938, 
+ 0.430060315238938, 0.430060315238938, 0.430060315238938, 0.430060315238938, 
+ 0.430060315238938, 0.430060315238938, 0.430060315238938, 0.430060315238938
+ ), INFIN = c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), CORREL = c(0.5, 
+ 0.5, 0.5, 0.5, 0.5, 0.5, -0.5, 0.5, -2.75206388903997e-16, -3.44007986129997e-16, 
+ -0.5, -4.12809583355996e-16, 0.499999999999999, -6.19214375033994e-16, 
+ 0.5, -0.5, -4.12809583355996e-16, -5.50412777807995e-16, 0.5, 
+ 0.5, 0.5, -1.37603194451999e-16, -0.5, 0.5, -2.75206388903997e-16, 
+ -0.5, 0.5, -1.37603194451999e-16, -6.88015972259993e-17, -0.5, 
+ -2.75206388903997e-16, 0.5, -0.5, -2.06404791677998e-16, 0.5, 
+ 0.5, 6.88015972259993e-17, 0, -0.5, 0.5, -6.88015972259993e-17, 
+ -0.5, 0.5, -0.5, 0.5), DELTA = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0
+ ), MAXPTS = 25000L, ABSEPS = 0.001, RELEPS = 0, error = NaN, 
+     value = NaN, inform = 0L), .Names = c("N", "NU", "LOWER", 
+ "UPPER", "INFIN", "CORREL", "DELTA", "MAXPTS", "ABSEPS", "RELEPS", 
+ "error", "value", "inform"))
> 
> RS <-
+ c(403L, 480L, 641015092L, 1848202935L, -2124158291L, -2116162620L, 
+ 1818211306L, -796165035L, -1592745489L, -483415562L, -77025504L, 
+ -1708531485L, 2015614337L, 1987179504L, 1442495118L, 792268281L, 
+ -362319989L, 798403514L, 744985276L, -385857025L, -2143185067L, 
+ -2123763604L, 1037640898L, -2059956451L, -1706850921L, 514864382L, 
+ -1711546696L, -338706277L, -1586794759L, -132177848L, 1947022198L, 
+ 1396713137L, 1307682531L, -1035135566L, -557613564L, -917080761L, 
+ -187518947L, -5106028L, 355507770L, -1769639803L, 277345631L, 
+ -1350443546L, 1186100336L, 1328736147L, 971714609L, 940736224L, 
+ -351651970L, -869421559L, -1443183045L, -895592438L, 512985068L, 
+ -1752126353L, 653214629L, 2037080924L, 819580722L, -188940755L, 
+ -1248906009L, -988531186L, 1246358504L, -1959914005L, -444776823L, 
+ 93416984L, 1671046470L, -737094751L, 923475507L, -1670573150L, 
+ 1363142292L, 1367254295L, 1064674061L, 203526692L, -730506742L, 
+ 1768186741L, 2072438479L, -54040106L, 1530643840L, 1545469379L, 
+ -1047274655L, -1790344752L, -929146578L, 988160217L, 1908423403L, 
+ 1908139866L, -1951620836L, 825419807L, 407434357L, 571175180L, 
+ -1050354462L, 357905853L, 2118664759L, 1473302750L, -1574953960L, 
+ 1687716731L, 830059289L, 980979624L, 1322192726L, 646707153L, 
+ -1536810109L, -60324334L, -815619676L, -1378118297L, 342137533L, 
+ 2121800244L, -1614339686L, 1127874021L, -1298114177L, 1086320902L, 
+ -404911472L, -2033165517L, 1667637457L, 610202240L, -612375714L, 
+ -1394552663L, -1349307557L, 1259094378L, -2052643828L, 1750860879L, 
+ -187222779L, -262015556L, 993968274L, 1228288909L, 1010942919L, 
+ 389707950L, 620378376L, 1577649931L, -1909325335L, -2390344L, 
+ 372404582L, 98136577L, 1088547987L, -721889470L, -1123440396L, 
+ -311384841L, 1731935469L, -183242876L, 1163420074L, -1161226475L, 
+ -1800960593L, 2010688950L, 1649464288L, 285670179L, -1872495551L, 
+ 249728816L, -1709367986L, 1208557497L, 2091882955L, 463919738L, 
+ -1663742084L, -2023678017L, 852953237L, -670894420L, 1470882946L, 
+ -1409563555L, -259257513L, 615437886L, 1163713144L, 933718363L, 
+ -2103899975L, -1780750072L, 1279171510L, -109061903L, 218861347L, 
+ 842111346L, 1183468484L, 1082697607L, -330715811L, 1795654484L, 
+ 1476237690L, -2071423803L, 875808287L, -1827581530L, 1367491376L, 
+ -1300303405L, 1167325169L, -1181200096L, -111890498L, -1160962103L, 
+ -1837420677L, 45192522L, 1886183852L, 1827503279L, -833690011L, 
+ -1287037412L, -95219598L, -593361171L, -62608601L, -1463236786L, 
+ -735514712L, -1059193941L, -2054260535L, -847052584L, 1957158278L, 
+ -158272415L, -2130332173L, -51186590L, -2085899564L, -575802409L, 
+ -853841459L, 1167301988L, -381429430L, -225222475L, 1258276879L, 
+ -313335530L, -1153729600L, -642721917L, -1628818015L, 276117904L, 
+ 1635992174L, 916804633L, 68094123L, -1364902758L, 1741683804L, 
+ -1146295969L, -2054907083L, 1265286860L, 1344646178L, 227339901L, 
+ 1398444919L, 684092574L, 1432542040L, 1360666299L, -1757052839L, 
+ 1498842088L, -616937962L, 1152706961L, 444549699L, -949761710L, 
+ 578780644L, -862362585L, 833589373L, 840321860L, -180485672L, 
+ -77392774L, 1715906096L, 1909965500L, 1289343828L, -1648824254L, 
+ 180038048L, -142870852L, 59178448L, 1076108642L, -483711000L, 
+ 1679486924L, -603249924L, -2138180110L, -1800352640L, 1949866228L, 
+ 299899528L, -1627927750L, 1867114320L, 1579244140L, 925202580L, 
+ -79527598L, -127183136L, -1460746548L, 1246824608L, -856848510L, 
+ 1497933544L, 1150267116L, -763654084L, 1936811634L, 204944496L, 
+ -134657500L, 2084542520L, -309653478L, 525739728L, 817298492L, 
+ 1309740180L, 490932354L, 856289216L, 2000386364L, -1149053296L, 
+ 1155498018L, -2127677432L, -1708438388L, 1378391068L, -1990972782L, 
+ 693629056L, 765580308L, 1863371048L, 99289658L, -341257008L, 
+ 949900396L, 841842356L, -1326807022L, 639191008L, 148198348L, 
+ 1457253856L, -1054367198L, -1883833048L, 102257932L, 1416220220L, 
+ 605913330L, 1216892656L, 991629124L, -2125930920L, 1549433402L, 
+ -1582915856L, -1950027588L, -1580588524L, 892478658L, 325637216L, 
+ 400973180L, -501153584L, -617483230L, 1942490792L, -1873501300L, 
+ -1073185348L, 377056306L, -417347776L, -1427909452L, -1015486648L, 
+ 1610431162L, -1914314992L, -973157140L, -564128684L, 850120466L, 
+ -480278624L, -434555956L, -996334752L, 1533891778L, -1993411032L, 
+ 1229952428L, -1688761540L, -1130449038L, 1459599152L, -332540380L, 
+ 870782008L, 1219178714L, -1356266096L, -1292801284L, -1685732460L, 
+ -1524381886L, -132346368L, -668895684L, 1288698448L, 612001442L, 
+ -620952760L, -319344116L, 157984988L, -1912149422L, -2012630912L, 
+ -1887089836L, -1111091736L, -775823622L, 1614588752L, 2116361260L, 
+ -830609036L, 1858582162L, -792424864L, 1982613132L, -1094686880L, 
+ -1498166558L, 882446504L, 102460620L, -1177968132L, 1414746034L, 
+ 92678640L, -1439926844L, 1724226392L, 381562746L, 289211824L, 
+ 1075324348L, 1606419540L, -761653310L, -587540704L, 1082387388L, 
+ -1655536048L, -821229086L, -736311960L, -1474683572L, -73616132L, 
+ 159378162L, -497990912L, -675148172L, 1318966280L, 195340730L, 
+ 16528336L, 475194860L, 1310052244L, -2015985966L, -582305184L, 
+ 1144814924L, -994553568L, -302756222L, 27157992L, 51995628L, 
+ 1116674620L, 500909682L, -1316848784L, 949899556L, -1881738696L, 
+ -556418662L, -86078000L, 1661341116L, 1430918292L, 182278274L, 
+ 1806048704L, 315494204L, 997703824L, 603310882L, 223906312L, 
+ -73062132L, 1177831836L, -1014101486L, -500732032L, 1207241492L, 
+ -1274969816L, -413595590L, 1742887888L, -1961988500L, 1086927412L, 
+ -1641077102L, 633855328L, -856995508L, 1927246816L, -248700510L, 
+ 525318568L, 1586427532L, 715701692L, -589647246L, 595708912L, 
+ -2060615100L, 2128179288L, -629423046L, -1817381776L, -153623492L, 
+ 634512660L, 559745730L, 1281177312L, -563131268L, 599345872L, 
+ -1708906078L, -452584408L, 1534535436L, -342831684L, 1837164338L, 
+ 1890266816L, 1054445620L, 1148866248L, -70563654L, -715414896L, 
+ 860946156L, -1200514476L, 1402941074L, -684954336L, 1116317132L, 
+ -1819848992L, -994736958L, -1403732952L, -2015169620L, -842944836L, 
+ 1209407474L, 733510192L, 915891108L, 809020856L, 1286688602L, 
+ -919972192L, -888651778L, 1926633163L, 438524109L, -1572701190L, 
+ 235874360L, 1760006625L, 1137293971L, 1528770740L, 1128096194L, 
+ -535410217L, -720898527L, -951685658L, -1292639004L, 1557699605L, 
+ -36201409L, -626756648L, 931948022L, 1093010947L, 1782374725L, 
+ -785382046L, 1155618512L, 2063380985L, -143388373L, 194689692L, 
+ 850628810L, 1920883135L, 1404267785L, -44388834L, -1847692340L, 
+ -702324323L, 722836519L, 1867891280L, 473042670L, 1189283739L, 
+ 2014475773L, -1801279254L, 1171364744L, -867510991L, 437669283L, 
+ -2050696796L, -344851374L, -881842265L, -1233521679L, 293668662L, 
+ 2089051092L, 813893989L, 704484719L, -517560408L, 998231302L, 
+ -1733749709L, -1037192555L, -840989966L, 622119744L, 344596649L, 
+ 736388379L, -1599397140L, 1932082938L, 325734191L, -443240647L, 
+ -1834520242L, 927092636L, -1675698931L, 909133687L, -1183726720L, 
+ -1425106466L, -1237459285L, -907075155L, 334411290L, -1003018024L, 
+ -473918143L, 768187891L, -275149804L, -1671032414L, 235177655L, 
+ 22335617L, 1414719622L, -1013482300L, -1246848523L, -1485122465L, 
+ 904693432L, -556791914L, -227484253L, -1578696347L, 1593534594L, 
+ -1118305808L, -1996105831L, -49671797L, 105245820L, -1524087126L, 
+ -1739448097L, -783939991L, -2143579010L, -1565295508L, 296688829L, 
+ 97826183L, 2121142640L, 736838670L, -1551402693L, -290765283L, 
+ 2108857034L, -151153176L, 538376913L, -1294006589L, 1415622084L, 
+ 2108467826L, -1698155705L, -1451947247L, 21739158L, 1179667444L, 
+ 1171699077L, -1024792369L, 1297433032L, -176160666L, 202118803L, 
+ 1908786677L, -1725081646L, -522804192L, -392088055L, -1557123269L, 
+ -950721908L, -2129250534L, 1855314319L, 818067161L, -732243602L, 
+ -1107827460L, -1727676819L, 474069527L, 1746923488L, 1548626878L, 
+ -452408949L, 1808062989L, -1826531782L, 1809743480L, 156939809L, 
+ -391608109L, -1388152844L, -1247715198L, 1121925399L, -541487135L, 
+ 1941702950L, -1821366236L, -163591467L, -903842177L, 420804376L, 
+ 2129622582L, -1830090557L, 1720263941L, -1852093278L, 141202064L, 
+ 1449426489L, -1242759445L, -1642184740L, -289816054L, 2115568895L, 
+ 1662521673L, -1272364322L, 908882060L, -1725650851L, -597965209L, 
+ -1566869616L, -1222206546L, 1890198107L, -664658371L, -1032011606L, 
+ 345071944L, -1002412687L, 599773923L, -795600412L, 1751993866L
+ )
> 
> 
> f <- function() {
+ 
+     error <- 0; value <- 0; inform <- 0
+     ret <- .C(C_mvtdst, N = as.integer(ret$N),  
+                         NU = as.integer(ret$NU),
+                         LOWER = as.double(ret$LOWER), 
+                         UPPER = as.double(ret$UPPER), 
+                         INFIN = as.integer(ret$INFIN),
+                         CORREL = as.double(ret$CORREL),
+                         DELTA = as.double(ret$DELTA),
+                         MAXPTS = as.integer(ret$MAXPTS),
+                         ABSEPS = as.double(ret$ABSEPS),
+                         RELEPS = as.double(ret$RELEPS),
+                         error = as.double(error),
+                         value = as.double(value),
+                         inform = as.integer(inform), RND = 1L)
+     ret
+ 
+ }
> 
> ### this special seed triggers the problem
> ### error and value are NaN (already in FORTRAN)
> .Random.seed <- RS
> # stopifnot(!is.na(f()$value)) ### .C does not work here
> 
> ### check tail with new quantile algorithm
> p <- .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 <- .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 - .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
> 
> ### package schwartz97
> qmvnorm(p = .5, tail = "lower", mean = c(6.75044368, 0.04996326), 
+         sigma = rbind(c(0.10260550, 0.02096418),
+                       c(0.02096418, 0.16049956)))$quantile
[1] 6.750444
> stint <- c(6.75044332319072, 6.75044368) ## with very narrow start interval
> qmvnorm(p = .5, tail = "lower", mean = c(6.75044368, 0.04996326),
+         sigma = rbind(c(0.10260550, 0.02096418),
+           c(0.02096418, 0.16049956)), interval=stint)$quantile
[1] 6.750444
> 
> ### qmvnorm and qmvt should stop if supplied covariance matrix
> ### is not positive semidefinite <Shiyang_Ma@URMC.Rochester.edu>
> 
> R2=matrix(c(0.7071068, 0.6924398, 0.7054602, 0.7054602, 0.6292745,
+             0.6924398, 0.7071068, 0.6909812, 0.6909712, 0.6128670,
+             0.7054602, 0.6909812, 0.7071068, 0.7071068, 0.6278091,
+             0.7054602, 0.6909712, 0.7071068, 0.7071068, 0.6278091,
+             0.6292745, 0.6128670, 0.6278091, 0.6278091, 0.7071068),ncol=5)
> 
> call <- try(qmvnorm(p=1-0.0001726701,
+                     mean=c(-0.8752332, -0.9487915, -0.9719237,
+                            -0.5855204, -0.9046457),
+                     sigma=R2,tail='lower.tail')$quantile,
+             silent=TRUE)
> 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.tail')$quantile,
+             silent=TRUE)
> inherits(call, "try-error")
[1] TRUE
> grepl("Covariance matrix not positive semidefinite", geterrmessage())
[1] TRUE
> 
> ### qmvnorm was wrong for the univariate setting; reported by Chen-Wei
> ### <cwliu@ntnu.edu.tw>
> 
> 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 
  4.166   0.051   4.212