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.
1 2 3 
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 tprobabilities 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 3dimensional problems one can also use the TVPACK
routines
described by Genz (2004), which only handles semiinfinite 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
semidefinite) 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 onesided
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 
http://www.sci.wsu.edu/math/faculty/genz/homepage
Genz, A. and Bretz, F. (1999), Numerical computation of multivariate tprobabilities 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 tprobabilities. Journal of Computational and Graphical Statistics, 11, 950–971.
Genz, A. (2004), Numerical computation of rectangular bivariate and trivariate normal and tprobabilities, 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. SpringerVerlag, 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 simulationbased multiple comparisons. Biometrics, 43, 913–928.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62  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=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)))

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.
Please suggest features or report bugs with the GitHub issue tracker.
All documentation is copyright its authors; we didn't write any of that.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.