Running 'bugfix-tests.R' [171s]
  Running 'dmvnorm-Ex.R' [156s]
  Running 'plmvnorm-Ex.R' [41s]
  Comparing 'plmvnorm-Ex.Rout' to 'plmvnorm-Ex.Rout.save' ... OK
  Running 'regtest-TVPACK.R' [8s]
  Comparing 'regtest-TVPACK.Rout' to 'regtest-TVPACK.Rout.save' ... OK
  Running 'regtest-scores.R' [187s]
  Running 'rmvnorm.R' [4s]
  Running 'slpmvnorm.R' [170s]
  Running 'test-getInt.R' [1s]
  Running 'test-noisy-root.R' [12s]
  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)
  > 
  > 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)))
  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 <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 
    10.03    0.23   10.25 
Running the tests in 'tests/dmvnorm-Ex.R' failed.
Complete output:
  > 
  > library("mvtnorm")
  > 
  > set.seed(29)
  > 
  > chk <- function(...) stopifnot(isTRUE(all.equal(..., check.attributes = FALSE)))
  > 
  > ### N samples with N different covariance matrices
  > 
  > thischeck <- expression({
  +   N <- 10
  +   J <- 5
  +   lt <- ltMatrices(matrix(runif(N * J * (J + c(-1, 1)[dg + 1L]) / 2) + 1, ncol = N), 
  +                    diag = dg)
  +   lt <- ltMatrices(lt, diag = dg, byrow = br)
  +   Z <- matrix(rnorm(N * J), ncol = N)
  +   Y <- solve(lt, Z)
  +   ll1 <- sum(dnorm(Mult(lt, Y), log = TRUE)) + sum(log(diagonals(lt)))
  + 
  +   S <- as.array(Tcrossprod(solve(lt)))
  +   ll2 <- sum(l2 <- sapply(1:N, function(i) mvtnorm:::dmvnorm(x = Y[,i], sigma = S[,,i], log = TRUE)))
  +   chk(ll1, ll2)
  + 
  +   l3 <- ldmvnorm(obs = Y, invchol = lt, logLik = FALSE)
  +   l4 <- ldmvnorm(obs = Y, chol = solve(lt), logLik = FALSE)
  + 
  +   chk(l2, l3)
  +   chk(l2, l4)
  + 
  +   ll1 <- sum(dnorm(Mult(lt[1,], Y), log = TRUE)) + N * sum(log(diagonals(lt[1,])))
  + 
  +   S <- as.array(Tcrossprod(solve(lt)))
  +   ll2 <- sum(l2 <- sapply(1:N, function(i) mvtnorm:::dmvnorm(x = Y[,i], sigma = S[,,1], log = TRUE)))
  +   chk(ll1, ll2)
  + 
  +   l3 <- ldmvnorm(obs = Y, invchol = lt[1,], logLik = FALSE)
  +   l4 <- ldmvnorm(obs = Y, chol = solve(lt[1,]), logLik = FALSE)
  + 
  +   chk(l2, l3)
  +   chk(l2, l4)
  + 
  +   ### check scores
  +   if (require("numDeriv", quietly = TRUE)) {
  + 
  +     f <- function(L) {
  +       L <- ltMatrices(L, diag = dg, byrow = br)
  +       ldmvnorm(obs = Y, invchol = L)
  +     }
  + 
  +     s0 <- grad(f, unclass(lt))
  +     s1 <- sldmvnorm(obs = Y, invchol = lt)
  + 
  +     chk(Lower_tri(ltMatrices(matrix(s0, ncol = N), diag = dg, byrow = br), diag = dg), 
  +         Lower_tri(s1$invchol, diag = dg))
  + 
  +     f <- function(L) {
  +       L <- ltMatrices(L, diag = dg, byrow = br)
  +       ldmvnorm(obs = Y, chol = L)
  +     }
  + 
  +     s0 <- grad(f, unclass(lt))
  +     s1 <- sldmvnorm(obs = Y, chol = lt)
  + 
  +     chk(Lower_tri(ltMatrices(matrix(s0, ncol = N), diag = dg, byrow = br), diag = dg), 
  +         Lower_tri(s1$chol, diag = dg))
  + 
  +     f <- function(x)
  +       ldmvnorm(obs = x, invchol = lt)
  + 
  +     s0 <- grad(f, Y)
  +     s1 <- sldmvnorm(obs = Y, invchol = lt)
  + 
  +     chk(matrix(s0, ncol = N), s1$obs)
  +   }
  + })
  > 
  > dg <- TRUE
  > br <- FALSE
  > eval(thischeck)
  > 
  > dg <- FALSE
  > br <- FALSE
  > eval(thischeck)
  > 
  > dg <- FALSE
  > br <- TRUE
  > eval(thischeck)
  > 
  > dg <- FALSE
  > br <- FALSE
  > eval(thischeck)
  > 
  > 
  > proc.time()
     user  system elapsed 
     2.85    0.12    2.92 
Running the tests in 'tests/regtest-scores.R' failed.
Complete output:
  > 
  > library("mvtnorm")
  > 
  > set.seed(29081975)
  > 
  > chk <- function(...) 
  +     stopifnot(all.equal(..., tol = 1e-5, check.attributes = FALSE))
  > 
  > EVAL <- function(...) {}
  > 
  > if (require("numDeriv", quietly = TRUE) && 
  +     require("qrng", quietly = TRUE))
  +   EVAL <- eval
  > 
  > N <- 10
  > M <- 10000
  > MM <- M / N
  > 
  > prb <- 1:3 / 4
  > 
  > ### chol 
  > thischeck <- expression({
  +   J <- cJ + dJ
  +   W <- NULL
  +   if (dJ > 1)
  +     W <- t(ghalton(M, d = dJ - 1))
  + 
  +   prm <- matrix(runif(J * (J - 1) / 2), ncol = 1)
  +   C <- ltMatrices(prm, byrow = BYROW)
  +   Z <- matrix(rnorm(J * N), ncol = N)
  +   Y <- Mult(C, Z)
  +   obs <- NULL
  +   if (cJ)
  +       obs <- Y[1:cJ,,drop = FALSE]
  +   lwr <- upr <- NULL
  +   if (dJ) {
  +     lwr <- t(apply(Y[cJ + (1:dJ),,drop = FALSE], 1, function(y) {
  +       qy <- quantile(y, prob = prb)
  +       c(-Inf, qy)[cut(y, breaks = c(-Inf, qy, Inf))]
  +     }))
  +     upr <- t(apply(Y[cJ + (1:dJ),,drop = FALSE], 1, function(y) {
  +       qy <- quantile(y, prob = prb)
  +       c(qy, Inf)[cut(y, breaks = c(-Inf, qy, Inf))]
  +     }))
  +   }
  + 
  +   ll <- function(prm) {
  +     C <- ltMatrices(prm, byrow = BYROW)
  +     -ldpmvnorm(obs = obs, lower = lwr, upper = upr, chol = C, M = MM, w = W)
  +   }
  + 
  +   sc <- function(prm) {
  +     C <- ltMatrices(prm, byrow = BYROW)
  +     ret <- sldpmvnorm(obs = obs, lower = lwr, upper = upr, chol = C,
  +                       M = MM, w = W)$chol
  +     -rowSums(Lower_tri(ret))
  +   }
  + 
  +   theta <- runif(J * (J - 1) / 2)
  +   print(ll(theta))
  +   chk(grad(ll, theta), sc(theta))
  + })
  > 
  > BYROW <- FALSE
  > cJ <- 4
  > dJ <- 4
  > EVAL(thischeck)
  [1] 110.5571
  > 
  > cJ <- 1
  > dJ <- 4
  > EVAL(thischeck)
  [1] 61.73499
  > 
  > cJ <- 4
  > dJ <- 1
  > EVAL(thischeck)
  [1] 101.5923
  > 
  > cJ <- 0
  > dJ <- 4
  > EVAL(thischeck)
  [1] 62.74771
  > 
  > cJ <- 4
  > dJ <- 0
  > EVAL(thischeck)
  [1] 68.67776
  > 
  > 
  > BYROW <- TRUE
  > cJ <- 4
  > dJ <- 4
  > EVAL(thischeck)
  [1] 126.3142
  > 
  > cJ <- 1
  > dJ <- 4
  > EVAL(thischeck)
  [1] 75.82315
  > 
  > cJ <- 4
  > dJ <- 1
  > EVAL(thischeck)
  [1] 69.16547
  > 
  > cJ <- 0
  > dJ <- 4
  > EVAL(thischeck)
  [1] 55.98251
  > 
  > cJ <- 4
  > dJ <- 0
  > EVAL(thischeck)
  [1] 55.17074
  > 
  > ### invchol
  > thischeck <- expression({
  +   J <- cJ + dJ
  +   W <- NULL
  +   if (dJ > 1)
  +     W <- t(ghalton(M, d = dJ - 1))
  + 
  +   prm <- matrix(runif(J * (J - 1) / 2), ncol = 1)
  +   C <- ltMatrices(prm, byrow = BYROW)
  +   Z <- matrix(rnorm(J * N), ncol = N)
  +   Y <- Mult(C, Z)
  +   obs <- NULL
  +   if (cJ)
  +       obs <- Y[1:cJ,,drop = FALSE]
  +   lwr <- upr <- NULL
  +   if (dJ) {
  +     lwr <- t(apply(Y[cJ + (1:dJ),,drop = FALSE], 1, function(y) {
  +       qy <- quantile(y, prob = prb)
  +       c(-Inf, qy)[cut(y, breaks = c(-Inf, qy, Inf))]
  +     }))
  +     upr <- t(apply(Y[cJ + (1:dJ),,drop = FALSE], 1, function(y) {
  +       qy <- quantile(y, prob = prb)
  +       c(qy, Inf)[cut(y, breaks = c(-Inf, qy, Inf))]
  +     }))
  +   }
  + 
  +   ll <- function(prm) {
  +     L <- ltMatrices(prm, byrow = BYROW)
  +     -ldpmvnorm(obs = obs, 
  +               lower = lwr, upper = upr, invchol = L, M = MM, w = W)
  +   }
  + 
  +   sc <- function(prm) {
  +     L <- ltMatrices(prm, byrow = BYROW)
  +     ret <- sldpmvnorm(obs = obs, 
  +                       lower = lwr, upper = upr, invchol = L,
  +                       M = MM, w = W)$invchol
  +     -rowSums(Lower_tri(ret))
  +   }
  + 
  +   theta <- runif(J * (J - 1) / 2)
  +   C <- ltMatrices(matrix(theta, ncol = 1), byrow = BYROW)
  +   theta <- Lower_tri(solve(C))
  +   print(ll(theta))
  +   chk(grad(ll, theta), sc(theta))
  + })
  > 
  > 
  > BYROW <- FALSE
  > cJ <- 4
  > dJ <- 4
  > EVAL(thischeck)
  [1] 119.9084
  > 
  > cJ <- 1
  > dJ <- 4
  > EVAL(thischeck)
  [1] 84.83088
  > 
  > cJ <- 4
  > dJ <- 1
  > EVAL(thischeck)
  [1] 79.17239
  > 
  > cJ <- 0
  > dJ <- 4
  > EVAL(thischeck)
  [1] 62.52789
  > 
  > cJ <- 4
  > dJ <- 0
  > EVAL(thischeck)
  [1] 70.16437
  > 
  > 
  > BYROW <- TRUE
  > cJ <- 4
  > dJ <- 4
  > EVAL(thischeck)
  [1] 125.3092
  > 
  > cJ <- 1
  > dJ <- 4
  > EVAL(thischeck)
  [1] 69.96545
  > 
  > cJ <- 4
  > dJ <- 1
  > EVAL(thischeck)
  [1] 74.51345
  > 
  > cJ <- 0
  > dJ <- 4
  > EVAL(thischeck)
  [1] 54.76586
  > 
  > cJ <- 4
  > dJ <- 0
  > EVAL(thischeck)
  [1] 56.95164
  > 
  > ### chol standardized
  > thischeck <- expression({
  +   J <- cJ + dJ
  +   W <- NULL
  +   if (dJ > 1)
  +     W <- t(ghalton(M, d = dJ - 1))
  + 
  +   prm <- matrix(runif(J * (J - 1) / 2), ncol = 1)
  +   C <- ltMatrices(prm)
  +   C <- ltMatrices(C, byrow = BYROW)
  +   Z <- matrix(rnorm(J * N), ncol = N)
  +   Y <- Mult(C, Z)
  +   obs <- NULL
  +   if (cJ)
  +       obs <- Y[1:cJ,,drop = FALSE]
  +   lwr <- upr <- NULL
  +   if (dJ) {
  +     lwr <- t(apply(Y[cJ + (1:dJ),,drop = FALSE], 1, function(y) {
  +       qy <- quantile(y, prob = prb)
  +       c(-Inf, qy)[cut(y, breaks = c(-Inf, qy, Inf))]
  +     }))
  +     upr <- t(apply(Y[cJ + (1:dJ),,drop = FALSE], 1, function(y) {
  +       qy <- quantile(y, prob = prb)
  +       c(qy, Inf)[cut(y, breaks = c(-Inf, qy, Inf))]
  +     }))
  +   }
  + 
  +   ll <- function(prm) {
  +     C <- ltMatrices(prm, byrow = BYROW)
  +     Cs <- standardize(chol = C)
  +     -ldpmvnorm(obs = obs, lower = lwr, upper = upr, chol = Cs, M = MM, w = W)
  +   }
  + 
  +   sc <- function(prm) {
  +     C <- ltMatrices(prm, byrow = BYROW)
  +     Cs <- standardize(chol = C)
  +     ret <- sldpmvnorm(obs = obs, lower = lwr, upper = upr, chol = Cs,
  +                       M = MM, w = W)$chol
  +     ret <- destandardize(chol = C, score_schol = ret)
  +     -rowSums(Lower_tri(ret))
  +   }
  + 
  +   theta <- runif(J * (J - 1) / 2)
  +   print(ll(theta))
  +   chk(grad(ll, theta), sc(theta))
  + })
  > 
  > BYROW <- FALSE
  > cJ <- 4
  > dJ <- 4
  > EVAL(thischeck)
  [1] 181.3173
  > 
  > cJ <- 1
  > dJ <- 4
  > EVAL(thischeck)
  [1] 64.84017
  > 
  > cJ <- 4
  > dJ <- 1
  > EVAL(thischeck)
  [1] 92.57033
  > 
  > cJ <- 0
  > dJ <- 4
  > EVAL(thischeck)
  [1] 57.32839
  > 
  > cJ <- 4
  > dJ <- 0
  > EVAL(thischeck)
  [1] 56.09638
  > 
  > 
  > BYROW <- TRUE
  > cJ <- 4
  > dJ <- 4
  > EVAL(thischeck)
  [1] 133.3763
  > 
  > cJ <- 1
  > dJ <- 4
  > EVAL(thischeck)
  [1] 83.61253
  > 
  > cJ <- 4
  > dJ <- 1
  > EVAL(thischeck)
  [1] 77.96918
  > 
  > cJ <- 0
  > dJ <- 4
  > EVAL(thischeck)
  [1] 56.77728
  > 
  > cJ <- 4
  > dJ <- 0
  > EVAL(thischeck)
  [1] 61.12252
  > 
  > ### invchol standardized
  > thischeck <- expression({
  +   J <- cJ + dJ
  +   W <- NULL
  +   if (dJ > 1)
  +     W <- t(ghalton(M, d = dJ - 1))
  + 
  +   prm <- matrix(runif(J * (J - 1) / 2), ncol = 1)
  +   C <- ltMatrices(prm, byrow = BYROW)
  +   Z <- matrix(rnorm(J * N), ncol = N)
  +   Y <- Mult(C, Z)
  +   obs <- NULL
  +   if (cJ)
  +       obs <- Y[1:cJ,,drop = FALSE]
  +   lwr <- upr <- NULL
  +   if (dJ) {
  +     lwr <- t(apply(Y[cJ + (1:dJ),,drop = FALSE], 1, function(y) {
  +       qy <- quantile(y, prob = prb)
  +       c(-Inf, qy)[cut(y, breaks = c(-Inf, qy, Inf))]
  +     }))
  +     upr <- t(apply(Y[cJ + (1:dJ),,drop = FALSE], 1, function(y) {
  +       qy <- quantile(y, prob = prb)
  +       c(qy, Inf)[cut(y, breaks = c(-Inf, qy, Inf))]
  +     }))
  +   }
  + 
  +   ll <- function(prm) {
  +     L <- ltMatrices(prm, byrow = BYROW)
  +     Ls <- standardize(invchol = L)
  +     -ldpmvnorm(obs = obs, 
  +               lower = lwr, upper = upr, invchol = Ls, M = MM, w = W)
  +   }
  + 
  +   sc <- function(prm) {
  +     L <- ltMatrices(prm, byrow = BYROW)
  +     Cs <- standardize(chol = solve(L))
  +     ret <- sldpmvnorm(obs = obs, 
  +                       lower = lwr, upper = upr, chol = Cs,
  +                       M = MM, w = W)$chol
  +     ret <- destandardize(invchol = L, score_schol = ret)
  +     -rowSums(Lower_tri(ret))
  +   }
  + 
  +   theta <- runif(J * (J - 1) / 2)
  +   C <- ltMatrices(matrix(theta, ncol = 1), byrow = BYROW)
  +   theta <- Lower_tri(solve(C))
  +   print(ll(theta))
  +   chk(grad(ll, theta), sc(theta))
  + })
  > 
  > 
  > BYROW <- FALSE
  > cJ <- 4
  > dJ <- 4
  > EVAL(thischeck)
  [1] 144.6606
  > 
  > cJ <- 1
  > dJ <- 4
  > EVAL(thischeck)
  [1] 99.0809
  > 
  > cJ <- 4
  > dJ <- 1
  > EVAL(thischeck)
  [1] 69.28283
  > 
  > cJ <- 0
  > dJ <- 4
  > EVAL(thischeck)
  [1] 52.48012
  > 
  > cJ <- 4
  > dJ <- 0
  > EVAL(thischeck)
  [1] 52.48746
  > 
  > 
  > BYROW <- TRUE
  > cJ <- 4
  > dJ <- 4
  > EVAL(thischeck)
  [1] 178.2047
  > 
  > cJ <- 1
  > dJ <- 4
  > EVAL(thischeck)
  [1] 98.60784
  > 
  > cJ <- 4
  > dJ <- 1
  > EVAL(thischeck)
  [1] 87.51641
  > 
  > cJ <- 0
  > dJ <- 4
  > EVAL(thischeck)
  [1] 49.60575
  > 
  > cJ <- 4
  > dJ <- 0
  > EVAL(thischeck)
  [1] 64.79248
  > 
  > proc.time()
     user  system elapsed 
    29.28    0.34   29.64 
Running the tests in 'tests/slpmvnorm.R' failed.
Complete output:
  > 
  > library("mvtnorm")
  > 
  > set.seed(290875)
  > 
  > if (require("numDeriv", quietly = TRUE)) {
  + 
  + 
  + chk <- function(...) stopifnot(isTRUE(all.equal(...)))
  + 
  + J <- 5
  + 
  + chks <- function(dg, tol = .Machine$double.eps^(1 / 4)) {
  + 
  +     prm <- runif(J * (J + c(-1, 1)[dg + 1L]) / 2)
  +     L <- ltMatrices(prm, diag = dg)
  +     a <- matrix(-2, nrow = J, ncol = 1)
  +     b <- matrix( 2, nrow = J, ncol = 1)
  +     M <- 10000L
  + 
  +     l <- function(x) {
  +         x <- ltMatrices(x, diag = dg)
  +         lpmvnorm(a, b, chol = x, M = M, seed = 29)
  +     }
  + 
  +     s <- function(x)
  +         slpmvnorm(a, b, chol = x, M = M, seed = 29)
  + 
  +     rl <- l(L)
  +     rs <- s(L)
  +     chk(rl, rs$logLik)
  +     if (dg) {
  +         chk(rs$chol, ltMatrices(grad(l, unclass(L)), diag = dg), tol = tol)
  +     } else {
  +         chk(c(Lower_tri(rs$chol)), grad(l, unclass(L)), tol = tol)
  +     }
  + 
  +     l <- function(x) {
  +         x <- ltMatrices(x, diag = dg)
  +         lpmvnorm(a, b, invchol = x, M = M, seed = 29)
  +     }
  + 
  +     s <- function(x)
  +         slpmvnorm(a, b, invchol = x, M = M, seed = 29)
  + 
  +     rl <- l(L)
  +     rs <- s(L)
  +     chk(rl, rs$logLik)
  +     if (dg) {
  +         chk(rs$invchol, ltMatrices(grad(l, unclass(L)), diag = dg), tol = tol)
  +     } else {
  +         chk(c(Lower_tri(rs$invchol)), grad(l, unclass(L)), tol = tol)
  +     }
  + 
  +     l <- function(x)
  +         lpmvnorm(a, b, mean = x, chol = L, M = M, seed = 29)
  +     s <- function(x)
  +         slpmvnorm(a, b, mean = x, chol = L, M = M, seed = 29)
  + 
  +     x <- numeric(J)
  +     rl <- l(x)
  +     rs <- s(x)
  +     chk(rl, rs$logLik)
  +     chk(grad(l, x), c(rs$mean), tol = tol)
  +     x <- 1:J
  +     rl <- l(x)
  +     rs <- s(x)
  +     chk(rl, rs$logLik)
  +     chk(grad(l, x), c(rs$mean), tol = tol)
  + }
  + 
  + chks(TRUE)
  + chks(FALSE)
  + 
  + ### check scores for conditional distributions
  + 
  + .cmvnorm <- function(invchol, which_given, given) {
  +     L <- invchol
  +     J <- dim(L)[2L]
  +     tmp <- matrix(0, ncol = ncol(given), nrow = J - length(which_given))
  +     center <- Mult(L, rbind(given, tmp))
  +     center <- center[-which_given,,drop = FALSE]
  +     L <- L[,-which_given]
  +     return(list(center = center, invchol = L))
  + }
  + 
  + J <- (cJ <- 4) + (dJ <- 4)
  + N <- 3
  + M <- 30
  + ltM <- function(x) ltMatrices(x, diag = FALSE, byrow = TRUE)
  + ltD <- function(x) ltMatrices(x, diag = TRUE, byrow = TRUE)
  + prm <- matrix(runif(J * (J - 1) / 2 * N), ncol = N)
  + L <- ltM(prm)
  + 
  + obs <- matrix(rnorm(J * N), ncol = N)
  + lwr <- -abs(obs)
  + upr <- abs(obs)
  + 
  + w <- matrix(runif((dJ - 1) * M), ncol = M)
  + 
  + ### without scaling, diag(L) == 1
  + ## score wrt L
  + j <- 1:cJ
  + ll <- function(x) {
  +     L <- ltM(x)
  +     cd <- .cmvnorm(invchol = L, which = j, given = obs[j,,drop = FALSE])
  +     lpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, 
  +             invchol = cd$invchol, w = w)
  + }
  + a <- ltM(matrix(grad(ll, unclass(L)), ncol = N))
  + diagonals(a) <- 0
  + 
  + cd <- .cmvnorm(invchol = L, which = j, given = obs[j,,drop = FALSE])
  + s <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, 
  +              invchol = cd$invchol, w = w)
  + 
  + chk(as.array(a[,-j]), as.array(s$invchol), check.attributes = FALSE)
  + 
  + ## score wrt obs
  + ll <- function(x) {
  +     cd <- .cmvnorm(invchol = L, which = j, given = x)
  +     lpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, 
  +             invchol = cd$invchol, w = w)
  + }
  + a <- matrix(grad(ll, obs[1:cJ,,drop = FALSE]), ncol = N)
  + 
  + cd <- .cmvnorm(invchol = L, which = j, given = obs[j,,drop = FALSE])
  + s <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, invchol = cd$invchol, w = w)
  + 
  + tmp0 <- solve(cd$invchol, s$mean, transpose = TRUE)
  + aL <- as.array(L)[-(1:cJ), 1:cJ,]
  + lst <- tmp0[rep(1:dJ, cJ),,drop = FALSE]
  + dobs <- -margin.table(aL * array(lst, dim = dim(aL)), 2:3)
  + 
  + chk(a, dobs, check.attributes = FALSE)
  + 
  + ## score wrt lower
  + ll <- function(x) {
  +     cd <- .cmvnorm(invchol = L, which = j, given = obs[1:cJ,,drop = FALSE])
  +     lpmvnorm(x, upr[-j,,drop = FALSE], center = cd$center, 
  +             invchol = cd$invchol, w = w)
  + }
  + a <- matrix(grad(ll, lwr[-(1:cJ),,drop = FALSE]), ncol = N)
  + 
  + cd <- .cmvnorm(invchol = L, which = j, given = obs[j,,drop = FALSE])
  + s <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, invchol = cd$invchol, w = w)
  + 
  + chk(a, s$lower, check.attributes = FALSE)
  + 
  + ## score wrt upper
  + ll <- function(x) {
  +     cd <- .cmvnorm(invchol = L, which = j, given = obs[1:cJ,,drop = FALSE])
  +     lpmvnorm(lwr[-j,,drop = FALSE], x, center = cd$center, 
  +             invchol = cd$invchol, w = w)
  + }
  + a <- matrix(grad(ll, upr[-(1:cJ),,drop = FALSE]), ncol = N)
  + 
  + cd <- .cmvnorm(invchol = L, which = j, given = obs[j,,drop = FALSE])
  + s <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, invchol = cd$invchol, w = w)
  + 
  + chk(a, s$upper, check.attributes = FALSE)
  + 
  + ### after scaling
  + LD <- invcholD(L)
  + 
  + ## score wrt LD (!)
  + # use center
  + ll <- function(x) {
  +     LD <- ltD(x)
  +     cd <- .cmvnorm(invchol = LD, which = j, given = obs[j,,drop = FALSE])
  +     lpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, 
  +             invchol = cd$invchol, w = w, logLik = TRUE)
  + }
  + 
  + a1 <- ltD(matrix(grad(ll, unclass(LD)), ncol = N))
  + # use mean
  + ll <- function(x) {
  +     LD <- ltD(x)
  +     cd <- cond_mvnorm(invchol = LD, which = j, given = obs[j,,drop = FALSE])
  +     lpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], mean = cd$mean, 
  +             invchol = cd$invchol, w = w, logLik = FALSE)
  + }
  + 
  + a2 <- ltMatrices(matrix(grad(function(...) sum(ll(...)), unclass(LD)), ncol = N), 
  +                 diag = TRUE, byrow = TRUE)
  + 
  + chk(a1, a2)
  + 
  + cd <- cond_mvnorm(invchol = LD, which = j, given = obs[j,,drop = FALSE])
  + s1 <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], mean = cd$mean, 
  +              invchol = cd$invchol, w = w)
  + 
  + cd <- .cmvnorm(invchol = LD, which = j, given = obs[j,,drop = FALSE])
  + s2 <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, 
  +              invchol = cd$invchol, w = w)
  + 
  + ### needs to be FALSE
  + # all.equal(s1, s2)
  + 
  + chk(a1[,-j], s2$invchol, check.attributes = FALSE)
  + 
  + # score wrt obs
  + ll <- function(x) {
  +     cd <- .cmvnorm(invchol = LD, which = j, given = x)
  +     lpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, 
  +             invchol = cd$invchol, w = w)
  + }
  + a <- matrix(grad(ll, obs[1:cJ,,drop = FALSE]), ncol = N)
  + 
  + cd <- .cmvnorm(invchol = LD, which = j, given = obs[j,,drop = FALSE])
  + s <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, invchol = cd$invchol, 
  +         w = w)
  + 
  + tmp0 <- solve(cd$invchol, s$mean, transpose = TRUE)
  + aL <- as.array(LD)[-(1:cJ), 1:cJ,]
  + lst <- tmp0[rep(1:dJ, cJ),,drop = FALSE]
  + dobs <- -margin.table(aL * array(lst, dim = dim(aL)), 2:3)
  + 
  + chk(a, dobs, check.attributes = FALSE)
  + 
  + ## score wrt lower
  + ll <- function(x) {
  +     cd <- .cmvnorm(invchol = LD, which = j, given = obs[1:cJ,,drop = FALSE])
  +     lpmvnorm(x, upr[-j,,drop = FALSE], center = cd$center, 
  +             invchol = cd$invchol, w = w)
  + }
  + a <- matrix(grad(ll, lwr[-(1:cJ),,drop = FALSE]), ncol = N)
  + 
  + cd <- .cmvnorm(invchol = LD, which = j, given = obs[j,,drop = FALSE])
  + s <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, invchol = cd$invchol, w = w)
  + 
  + chk(a, s$lower, check.attributes = FALSE)
  + 
  + ## score wrt upper
  + ll <- function(x) {
  +     cd <- .cmvnorm(invchol = LD, which = j, given = obs[1:cJ,,drop = FALSE])
  +     lpmvnorm(lwr[-j,,drop = FALSE], x, center = cd$center, 
  +             invchol = cd$invchol, w = w)
  + }
  + a <- matrix(grad(ll, upr[-(1:cJ),, drop = FALSE]), ncol = N)
  + 
  + cd <- .cmvnorm(invchol = LD, which = j, given = obs[j,,drop = FALSE])
  + s <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, invchol = cd$invchol, w = w)
  + 
  + chk(a, s$upper, check.attributes = FALSE)
  + 
  + ### one-dimensional conditional distribution
  + 
  + J <- (cJ <- 4) + (dJ <- 1)
  + prm <- matrix(runif(J * (J - 1) / 2 * N), ncol = N)
  + L <- ltM(prm)
  + 
  + obs <- matrix(rnorm(J * N), ncol = N)
  + lwr <- -abs(obs)
  + upr <- abs(obs)
  + 
  + w <- matrix(runif((dJ - 1) * M), ncol = M)
  + 
  + ### without scaling, diag(L) == 1
  + ## score wrt L not needed (a constant)
  + j <- 1:cJ
  + 
  + ## score wrt obs
  + ll <- function(x) {
  +     cd <- .cmvnorm(invchol = L, which = j, given = x)
  +     lpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, 
  +             invchol = cd$invchol, w = w)
  + }
  + a <- matrix(grad(ll, obs[1:cJ,,drop = FALSE]), ncol = N)
  + 
  + cd <- .cmvnorm(invchol = L, which = j, given = obs[j,,drop = FALSE])
  + s <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, invchol = cd$invchol, w = w)
  + 
  + tmp0 <- solve(cd$invchol, s$mean, transpose = TRUE)
  + aL <- as.array(L)[-(1:cJ), 1:cJ,]
  + lst <- tmp0[rep(1:dJ, cJ),,drop = FALSE]
  + dobs <- -aL * array(lst, dim = dim(aL))
  + 
  + chk(a, dobs, check.attributes = FALSE)
  + 
  + ## score wrt lower
  + ll <- function(x) {
  +     cd <- .cmvnorm(invchol = L, which = j, given = obs[1:cJ,,drop = FALSE])
  +     lpmvnorm(x, upr[-j,,drop = FALSE], center = cd$center, 
  +             invchol = cd$invchol, w = w)
  + }
  + a <- matrix(grad(ll, lwr[-(1:cJ),,drop = FALSE]), ncol = N)
  + 
  + cd <- .cmvnorm(invchol = L, which = j, given = obs[j,,drop = FALSE])
  + s <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, invchol = cd$invchol, w = w)
  + 
  + chk(a, s$lower, check.attributes = FALSE)
  + 
  + ## score wrt upper
  + ll <- function(x) {
  +     cd <- .cmvnorm(invchol = L, which = j, given = obs[1:cJ,,drop = FALSE])
  +     lpmvnorm(lwr[-j,,drop = FALSE], x, center = cd$center, 
  +             invchol = cd$invchol, w = w)
  + }
  + a <- matrix(grad(ll, upr[-(1:cJ),,drop = FALSE]), ncol = N)
  + 
  + cd <- .cmvnorm(invchol = L, which = j, given = obs[j,,drop = FALSE])
  + s <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, invchol = cd$invchol, w = w)
  + 
  + chk(a, s$upper, check.attributes = FALSE)
  + 
  + ### after scaling
  + LD <- invcholD(L)
  + 
  + ## score wrt LD (!)
  + # use center
  + ll <- function(x) {
  +     LD <- ltD(x)
  +     cd <- .cmvnorm(invchol = LD, which = j, given = obs[j,,drop = FALSE])
  +     lpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, 
  +             invchol = cd$invchol, w = w, logLik = TRUE)
  + }
  + 
  + a1 <- ltD(matrix(grad(ll, unclass(LD)), ncol = N))
  + # use mean
  + ll <- function(x) {
  +     LD <- ltD(x)
  +     cd <- cond_mvnorm(invchol = LD, which = j, given = obs[j,,drop = FALSE])
  +     lpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], mean = cd$mean, 
  +             invchol = cd$invchol, w = w, logLik = FALSE)
  + }
  + 
  + a2 <- ltMatrices(matrix(grad(function(...) sum(ll(...)), unclass(LD)), ncol = N), 
  +                 diag = TRUE, byrow = TRUE)
  + 
  + chk(a1, a2)
  + 
  + cd <- cond_mvnorm(invchol = LD, which = j, given = obs[j,,drop = FALSE])
  + s1 <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], mean = cd$mean, 
  +              invchol = cd$invchol, w = w)
  + 
  + cd <- .cmvnorm(invchol = LD, which = j, given = obs[j,,drop = FALSE])
  + s2 <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, 
  +              invchol = cd$invchol, w = w)
  + 
  + ### needs to be FALSE
  + # all.equal(s1, s2)
  + 
  + chk(a1[,-j], s2$invchol, check.attributes = FALSE)
  + 
  + # score wrt obs
  + ll <- function(x) {
  +     cd <- .cmvnorm(invchol = LD, which = j, given = x)
  +     lpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, 
  +             invchol = cd$invchol, w = w)
  + }
  + a <- matrix(grad(ll, obs[1:cJ,,drop = FALSE]), ncol = N)
  + 
  + cd <- .cmvnorm(invchol = LD, which = j, given = obs[j,,drop = FALSE])
  + s <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, invchol = cd$invchol, 
  +         w = w)
  + 
  + tmp0 <- solve(cd$invchol, s$mean, transpose = TRUE)
  + aL <- as.array(LD)[-(1:cJ), 1:cJ,]
  + lst <- tmp0[rep(1:dJ, cJ),,drop = FALSE]
  + dobs <- -(aL * array(lst, dim = dim(aL)))
  + 
  + chk(a, dobs, check.attributes = FALSE)
  + 
  + ## score wrt lower
  + ll <- function(x) {
  +     cd <- .cmvnorm(invchol = LD, which = j, given = obs[1:cJ,,drop = FALSE])
  +     lpmvnorm(x, upr[-j,,drop = FALSE], center = cd$center, 
  +             invchol = cd$invchol, w = w)
  + }
  + a <- matrix(grad(ll, lwr[-(1:cJ),,drop = FALSE]), ncol = N)
  + 
  + cd <- .cmvnorm(invchol = LD, which = j, given = obs[j,,drop = FALSE])
  + s <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, invchol = cd$invchol, w = w)
  + 
  + chk(a, s$lower, check.attributes = FALSE)
  + 
  + ## score wrt upper
  + ll <- function(x) {
  +     cd <- .cmvnorm(invchol = LD, which = j, given = obs[1:cJ,,drop = FALSE])
  +     lpmvnorm(lwr[-j,,drop = FALSE], x, center = cd$center, 
  +             invchol = cd$invchol, w = w)
  + }
  + a <- matrix(grad(ll, upr[-(1:cJ),,drop = FALSE]), ncol = N)
  + 
  + cd <- .cmvnorm(invchol = LD, which = j, given = obs[j,,drop = FALSE])
  + s <- slpmvnorm(lwr[-j,,drop = FALSE], upr[-j,,drop = FALSE], center = cd$center, invchol = cd$invchol, w = w)
  + 
  + chk(a, s$upper, check.attributes = FALSE)
  + 
  + ### check scores for extremely small likelihoods
  + ### use independence model as ground truth
  + J <- 10
  + 
  + yl <- matrix(-as.double(rep(1, J) / 2000), ncol = 1)
  + yr <- abs(yl)
  + 
  + w <- matrix(.5, nrow = J - 1, ncol = 1)
  + 
  + ## L = diag(J)
  + L <- ltMatrices(rep(1, J * (J + 1) / 2), diag = TRUE)
  + L[] <- 0
  + diagonals(L) <- 1
  + 
  + lpmvnorm(lower = yl, upper = yr, invchol = L, w = w)
  + 
  + s <- slpmvnorm(lower = yl, upper = yr, invchol = L, w = w)[c("logLik", "invchol")]
  + 
  + chk(c((dnorm(yr) * yr - dnorm(yl) * yl ) / (pnorm(yr) - pnorm(yl))),
  +     c(diagonals(s$invchol)))
  + }
  > 
  > proc.time()
     user  system elapsed 
    12.79    0.12   12.93