pmvt {mvtnorm} R Documentation

## Multivariate t Distribution

### Description

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.

### Usage

```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)
```

### Arguments

 `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.

### Details

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.

### Value

The evaluated distribution function is returned with attributes

 `error` estimated absolute error and `msg` status messages.

### Author(s)

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>

### References

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.

`qmvt`

### Examples

```
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))

```

[Package mvtnorm version 0.7-1 Index]