Multivariate t Distribution
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 function 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 (for default algorithm =
GenzBretz()) is based on randomized quasi Monte Carlo methods and
described in Genz and Bretz (1999, 2002).
Because of the randomization, the result for this algorithm (slightly)
depends on .Random.seed.
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),
mode (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 message (a |
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
}
### uniroot for this simple problem
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=300,
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)))Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.