Running ‘bugfix-tests.R’ [12s/14s]
  Running ‘regtest-TVPACK.R’ [9s/10s]
  Comparing ‘regtest-TVPACK.Rout’ to ‘regtest-TVPACK.Rout.save’ ... OK
  Running ‘rmvnorm.R’ [5s/6s]
  Running ‘test-getInt.R’ [1s/2s]
  Running ‘test-noisy-root.R’ [13s/14s]
  Comparing ‘test-noisy-root.Rout’ to ‘test-noisy-root.Rout.save’ ... OK
Running the tests in ‘tests/bugfix-tests.R’ failed.
Complete output:
  > invisible(options(echo = TRUE))
  > library("mvtnorm")
  > set.seed(290875)
  > 
  > # 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(all.equal(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(all.equal(a, pnorm(2, sd=sqrt(1.5))))
  > a <- pmvnorm(lower=-Inf,upper=2,mean=0,sigma=matrix(.5))
  > attributes(a) <- NULL
  > stopifnot(all.equal(a, pnorm(2, sd=sqrt(.5))))
  > a <- pmvnorm(lower=-Inf,upper=2,mean=0,sigma=.5)
  > attributes(a) <- NULL
  > stopifnot(all.equal(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(all.equal(dmvnorm(x=0, mean=30, log=TRUE),
  + 		    dnorm  (0,	      30, log=TRUE)))
  > 
  > stopifnot(
  +     all.equal(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))
  +     ,
  +     all.equal(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(all.equal(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(pmvnorm(upper=c(-Inf,1)) == 0)
  > stopifnot(pmvnorm(lower=c(Inf,1)) == 0)
  > stopifnot(pmvnorm(lower=c(-2,0),upper=c(-1,1),corr=matrix(rep(1,4),2,2)) == 0)
  > 
  > # bugged Fritz (long time ago)
  > stopifnot(all.equal(pmvnorm(-Inf, c(Inf, 0), 0, diag(2)),
  + 		    pmvnorm(-Inf, c(Inf, 0), 0)))
  > 
  > # this is a bug in `mvtdst' nobody was able to fix yet :-(
  > stopifnot(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(all.equal(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(all.equal(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(all.equal(a, pnorm(2, sd=sqrt(1.5))))
  > a <- pmvnormM(lower=-Inf,upper=2,mean=0,sigma=matrix(.5))
  > attributes(a) <- NULL
  > stopifnot(all.equal(a, pnorm(2, sd=sqrt(.5))))
  > a <- pmvnormM(lower=-Inf,upper=2,mean=0,sigma=.5)
  > attributes(a) <- NULL
  > stopifnot(all.equal(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(pmvnormM(upper=c(-Inf,1)) == 0)
  > stopifnot(pmvnormM(lower=c(Inf,1)) == 0)
  > 
  > # bugged Fritz (long time ago)
  > stopifnot(all.equal(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(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(all.equal(p, ptmp))
  > 
  > ### was == 1; spotted by Alex Lenkoski <lenkoski@stat.washington.edu>
  > stopifnot(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(all.equal(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
  
  $f.quantile
  [1] 0
  
  attr(,"message")
  [1] "Normal Completion"
  > 
  > ### 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(all.equal(ld1, ld2))
  + }
  Error: ld1 and ld2 are not equal:
    Mean relative difference: 4.895976e-07
  Execution halted