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, algorithm = GenzBretz(), type = c("Kshirsagar", "shifted"), ...)
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, for
|
df |
degree of freedom as integer. Normal probabilities are computed for |
corr |
the correlation matrix of dimension n. |
sigma |
the scale matrix of dimension n. Either |
algorithm |
an object of class |
type |
type of the noncentral multivariate t distribution
to be computed. |
... |
additional parameters (currently given to |
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 based on randomized quasi Monte Carlo methods and described in Genz and Bretz (1999, 2002).
For 2- and 3-dimensional problems one can also use the TVPACK routines
described by Genz (2004), which only handles semi-infinite integration
regions (and for type = "Kshirsagar"
only central problems).
For type = "Kshirsagar"
and a given correlation matrix
corr
, for short A, say, (which has to be positive
semi-definite) and degrees of freedom ν the following values are
numerically evaluated
I = 2^{1-ν/2} / Γ(ν/2) \int_0^∞ s^{ν-1} \exp(-s^2/2) Φ(s \cdot lower/√{ν} - δ, s \cdot upper/√{ν} - δ) \, ds
where
Φ(a,b) = (det(A)(2π)^m)^{-1/2} \int_a^b \exp(-x^\prime Ax/2) \, dx
is the multivariate normal distribution and m is the number of rows of A.
For type = "shifted"
, a positive definite symmetric matrix
S (which might be the correlation or the scale matrix),
non-centrality vector δ and degrees of freedom ν the
following integral is evaluated:
c\int_{lower_1}^{upper_1}...\int_{lower_m}^{upper_m} (1+(x-δ)'S^{-1}(x-δ)/ν)^{-(ν+m)/2}\, dx_1 ... dx_m,
where
c = Γ((ν+m)/2)/((π ν)^{m/2}Γ(ν/2)|S|^{1/2}),
and m is the number of rows of S.
Note that both -Inf
and +Inf
may be specified in
the lower and upper integral limits in order to compute one-sided
probabilities.
Univariate problems are passed to pt
.
If df = 0
, normal probabilities are returned.
The evaluated distribution function is returned with attributes
error |
estimated absolute error and |
msg |
status messages. |
http://www.sci.wsu.edu/math/faculty/genz/homepage
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.
Genz, A. (2004), Numerical computation of rectangular bivariate and trivariate normal and t-probabilities, Statistics and Computing, 14, 251–260.
Genz, A. and Bretz, F. (2009), Computation of Multivariate Normal and t Probabilities. Lecture Notes in Statistics, Vol. 195. Springer-Verlag, Heidelberg.
S. Kotz and S. Nadarajah (2004), Multivariate t Distributions and Their Applications. Cambridge University Press. Cambridge.
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) ### scale 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 scale 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)))