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. > 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 > > 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 > 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 > 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 > # 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 > 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 > 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 > 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 > 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))) 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 > > # 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 > 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 > 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 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 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 > 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 > stopifnot(is.finite(qmvt(.95, df = 0, corr = matrix(1))$quantile)) > > ### spotted by 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 > ### 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 > > 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 > ### > > 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 11.48 0.18 11.62