Nothing
tests/bugfix-tests.R
In mvtnorm: Multivariate Normal and t Distributions
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)
# 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)
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 <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
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 <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)
# 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)))
# 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
pmvt(lower, upper, sigma = sigma, delta=mn, df=df, type = "shifted")
ind2 <- apply(sims2, 1, comp, lower=lower, upper=upper)
mean(ind2)
pmvt(lower, upper, sigma = sigma, delta=mn, df=df, type = "Kshirsagar")
## 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
pmvt(lower, upper, corr = corrMat, delta=mn, df=df, type = "shifted")
ind2 <- apply(sims2, 1, comp, lower=lower, upper=upper)
mean(ind2)
pmvt(lower, upper, corr = corrMat, delta=mn, df=df, type = "Kshirsagar")
### 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
### 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
(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 <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
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 <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")
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
### <cwliu@ntnu.edu.tw>
all.equal(qnorm(p = 0.2397501, mean = 1, sd = sqrt(2)),
qmvnorm(p=0.2397501 , mean = 1, sigma = 2)$quantile)
Try the mvtnorm package in your browser
Any scripts or data that you put into this service are public.
mvtnorm documentation built on Sept. 11, 2024, 6:07 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
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)
# 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)
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 <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
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 <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)
# 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)))
# 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
pmvt(lower, upper, sigma = sigma, delta=mn, df=df, type = "shifted")
ind2 <- apply(sims2, 1, comp, lower=lower, upper=upper)
mean(ind2)
pmvt(lower, upper, sigma = sigma, delta=mn, df=df, type = "Kshirsagar")
## 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
pmvt(lower, upper, corr = corrMat, delta=mn, df=df, type = "shifted")
ind2 <- apply(sims2, 1, comp, lower=lower, upper=upper)
mean(ind2)
pmvt(lower, upper, corr = corrMat, delta=mn, df=df, type = "Kshirsagar")
### 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
### 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
(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 <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
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 <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")
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
### <cwliu@ntnu.edu.tw>
all.equal(qnorm(p = 0.2397501, mean = 1, sd = sqrt(2)),
qmvnorm(p=0.2397501 , mean = 1, sigma = 2)$quantile)
Try the mvtnorm package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.