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) # 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) dmvnorm(x=c(0,0), mean=c(25,25),log=TRUE) dmvnorm(x=c(0,0), mean=c(30,30),log=TRUE) 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 pmvt(lower=c(-Inf,-Inf), upper=c(2,2), delta=c(0, 0), df=25, corr=diag(2)) pmvt(lower=c(-Inf,-Inf), upper=c(2,2), delta=c(0, 0), df=250, corr=diag(2)) pmvt(lower=c(-Inf,-Inf), upper=c(2,2), delta=c(0, 0), df=1340, corr=diag(2)) pmvt(lower=c(-Inf,-Inf), upper=c(2,2), delta=c(0, 0), df=2500, corr=diag(2)) pmvt(lower=c(-100,-100), upper=c(2,2), delta=c(0, 0), df=2500, corr=diag(2)) # df = 0 pmvt(lower=c(-Inf,-Inf), upper=c(2,2), delta=c(0, 0), df=0, corr=diag(2)) pmvt(lower=-Inf, upper = 2, delta=0, df=0, corr=1) pnorm(2) # larger dimensions pnorm(2)^2 pmvnorm(lower=rep(-Inf, 2), upper=rep(2,2), sigma = diag(2)) pnorm(2)^90 pmvnorm(lower=rep(-Inf, 90), upper=rep(2,90), sigma = diag(90)) pnorm(2)^199 pmvnorm(lower=rep(-Inf, 199), upper=rep(2,199), sigma = diag(199)) # 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 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 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) # 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))) # 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(xlower) & all(x 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 ### 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 (Sig <- tcrossprod(L)) 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 qmvnorm(p, sigma = diag(d), tail = "upper")$quantile qmvnorm(p, sigma = diag(d), tail = "both")$quantile p <- 1 - .95 d <- 4 qmvnorm(p, sigma = diag(d), tail = "lower")$quantile qmvnorm(p, sigma = diag(d), tail = "upper")$quantile ### 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 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 ### 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") grepl("Covariance matrix not positive semidefinite", geterrmessage()) 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") grepl("Covariance matrix not positive semidefinite", geterrmessage()) ### 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)