| pmvt {mvtnorm} | R Documentation |
Computes the the distribution function of the multivariate t distribution for arbitrary limits, degrees of freedom and correlation matrices based on algorithms by Genz and Bretz.
pmvt(lower=-Inf, upper=Inf, delta=rep(0, length(lower)),
df=1, corr=NULL, sigma=NULL, maxpts = 25000, abseps = 0.001,
releps = 0)
rmvt(n, sigma=diag(2), df=1)
lower |
the vector of lower limits of length n. |
upper |
the vector of upper limits of length n. |
delta |
the vector of noncentrality parameters of length n. |
df |
degree of freedom as integer. |
corr |
the correlation matrix of dimension n. |
sigma |
the covariance matrix of dimension n. Either corr or
sigma can be specified. If sigma is given, the
problem is standardized. If neither corr nor
sigma is given, the identity matrix is used
for sigma. |
maxpts |
maximum number of function values as integer. |
abseps |
absolute error tolerance as double. |
releps |
relative error tolerance as double. |
n |
number of observations. |
This program involves the computation of central and noncentral multivariate t-probabilities with arbitrary correlation matrices. It involves both the computation of singular and nonsingular probabilities. The methodology is described in Genz and Bretz (1999, 2002).
For a given correlation matrix corr, for short A, say,
(which has to be positive semi-definite) and
degrees of freedom df the following
values are numerically evaluated
I = K int s^{df-1} exp(-s^2/2) Phi(s cdot lower/sqrt{df}-delta, s cdot upper/sqrt{df}-delta) ds
where Phi(a,b) = K^prime int_a^b exp(-x^prime Ax/2) dx is the multivariate normal distribution, K^prime = 1/sqrt{det(A)(2π)^m} and K = 2^{1-df/2} / Gamma(df/2) are constants and the (single) integral of I goes from 0 to +Inf.
Note that both -Inf and +Inf may be specified in
the lower and upper integral limits in order to compute one-sided
probabilities. Randomized quasi-Monte Carlo methods are used for the
computations.
Univariate problems are passed to pt.
Further information can be obtained from the quoted articles, which can be downloaded (together with additional material and additional codes) from the websites http://www.bioinf.uni-hannover.de/~bretz/ and http://www.sci.wsu.edu/math/faculty/genz/homepage.
rmvt is a wrapper to rmvnorm for random number
generation.
If df = 0, normal probabilities are returned.
The evaluated distribution function is returned with attributes
error |
estimated absolute error and |
msg |
status messages. |
Fortran Code by Alan Genz <AlanGenz@wsu.edu> and Frank Bretz <frank.bretz@pharma.novartis.com>, R port by Torsten Hothorn <Torsten.Hothorn@rzmail.uni-erlangen.de>
Genz, A. and Bretz, F. (1999), Numerical computation of multivariate t-probabilities with application to power calculation of multiple contrasts. Journal of Statistical Computation and Simulation, 63, 361–378.
Genz, A. and Bretz, F. (2002), Methods for the computation of multivariate t-probabilities. Journal of Computational and Graphical Statistics, 11, 950–971.
Edwards D. and Berry, Jack J. (1987), The efficiency of simulation-based multiple comparisons. Biometrics, 43, 913–928.
n <- 5
lower <- -1
upper <- 3
df <- 4
corr <- diag(5)
corr[lower.tri(corr)] <- 0.5
delta <- rep(0, 5)
prob <- pmvt(lower=lower, upper=upper, delta=delta, df=df, corr=corr)
print(prob)
pmvt(lower=-Inf, upper=3, df = 3, sigma = 1) == pt(3, 3)
# Example from R News paper (original by Edwards and Berry, 1987)
n <- c(26, 24, 20, 33, 32)
V <- diag(1/n)
df <- 130
C <- c(1,1,1,0,0,-1,0,0,1,0,0,-1,0,0,1,0,0,0,-1,-1,0,0,-1,0,0)
C <- matrix(C, ncol=5)
### covariance matrix
cv <- C %*% V %*% t(C)
### correlation matrix
dv <- t(1/sqrt(diag(cv)))
cr <- cv * (t(dv) %*% dv)
delta <- rep(0,5)
myfct <- function(q, alpha) {
lower <- rep(-q, ncol(cv))
upper <- rep(q, ncol(cv))
pmvt(lower=lower, upper=upper, delta=delta, df=df,
corr=cr, abseps=0.0001) - alpha
}
round(uniroot(myfct, lower=1, upper=5, alpha=0.95)$root, 3)
# compare pmvt and pmvnorm for large df:
a <- pmvnorm(lower=-Inf, upper=1, mean=rep(0, 5), corr=diag(5))
b <- pmvt(lower=-Inf, upper=1, delta=rep(0, 5), df=rep(300,5),
corr=diag(5))
a
b
stopifnot(round(a, 2) == round(b, 2))
# correlation and covariance matrix
a <- pmvt(lower=-Inf, upper=2, delta=rep(0,5), df=3,
sigma = diag(5)*2)
b <- pmvt(lower=-Inf, upper=2/sqrt(2), delta=rep(0,5),
df=3, corr=diag(5))
attributes(a) <- NULL
attributes(b) <- NULL
a
b
stopifnot(all.equal(round(a,3) , round(b, 3)))
a <- pmvt(0, 1,df=10)
attributes(a) <- NULL
b <- pt(1, df=10) - pt(0, df=10)
stopifnot(all.equal(round(a,10) , round(b, 10)))
rmvt(10, sigma=diag(10))