R version 3.1.0 (2014-04-10) -- "Spring Dance" Copyright (C) 2014 The R Foundation for Statistical Computing Platform: x86_64-unknown-linux-gnu (64-bit) R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > invisible(options(echo = TRUE)) > library("mvtnorm") > set.seed(290875) > > # 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(all.equal(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(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 > 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 > # 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 > 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 > 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 > 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 > 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 > 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 > 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 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)$quantile > q2 <- qmvnorm((k*(k-1))/(m*(m-1))*alpha, tail="upper", sigma=var, + interval = c(0, 5))$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.233641 $f.quantile [1] 1.652806e-07 attr(,"error") [1] 0 attr(,"msg") [1] "Normal Completion" $iter [1] 13 $init.it [1] NA $estim.prec [1] 6.103516e-05 > > ### 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(all.equal(ld1, ld2)) + } > > ### --- 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(identical(fx, c(Inf, rep(0, 1+10)))) > ## gave NaN for a long time, then error, then NaN, now we have converged ;-) > > > proc.time() user system elapsed 6.788 0.036 6.828