The multivariate Student's t distribution

The probability density function, the distribution function and random number generation for a d-dimensional Student's t random variable.

Usage

dmt(x, mean = rep(0, d), S, df=Inf, log = FALSE) 
pmt(x, mean = rep(0, d), S, df=Inf, ...) 
rmt(n = 1, mean = rep(0, d), S, df=Inf, sqrt=NULL) 
sadmvt(df, lower, upper, mean, S, maxpts = 2000*d, abseps = 1e-06, releps = 0) 
biv.nt.prob(df, lower, upper, mean, S)
ptriv.nt(df, x, mean, S)

Arguments

x: either a vector of length d or (for dmt and pmt) a matrix with d columns representing the coordinates of the point(s) where the density must be evaluated; see also ‘Details’.
mean: either a vector of length d, representing the location parameter (equal to the mean vector when df>1), or (for dmt and pmt) a matrix whose rows represent different mean vectors; in the matrix case, its dimensions must match those of x.
S: a symmetric positive definite matrix with dimensions (d,d) representing the scale matrix of the distribution, such that S*df/(df-2) is the variance-covariance matrix when df>2; a vector of length 1 is also allowed (in this case, d=1 is set).
df: the degrees of freedom. For rmt, it must be a positive real value or Inf. For all other functions, it must be a positive integer or Inf. A value df=Inf is translated to a call to a suitable function for the the multivariate normal distribution. See ‘Details’ for its effect for the evaluation of distribution functions and other probabilities.
log: a logical value(default value is FALSE); if TRUE, the logarithm of the density is computed.
sqrt: if not NULL (default value is NULL), a square root of the intended scale matrix S; see ‘Details’ for a full description.
...: arguments passed to sadmvt, among maxpts, absrel, releps.
n: the number of random vectors to be generated
lower: a numeric vector of lower integration limits of the density function; must be of maximal length 20; +Inf and -Inf entries are allowed.
upper: a numeric vector of upper integration limits of the density function; must be of maximal length 20; +Inf and -Inf entries are allowed
maxpts: the maximum number of function evaluations (default value: 2000*d)
abseps: absolute error tolerance (default value: 1e-6).
releps: relative error tolerance (default value: 0).

Details

The dimension d cannot exceed 20 for pmt and sadmvt. If this threshold is exceeded, NA is returned.

The functions sadmvt, ptriv.mt and biv.nt.prob are interfaces to Fortran 77 routines by Alan Genz, dowloaded from his web page; they makes use of some auxiliary functions whose authors are indicated in the Fortran code itself. The routine sadmvt uses an adaptive integration method. If df=3, a call to pmt activates a call to ptriv.nt which is specific for the trivariate case, and uses Genz's Fortran code tvpack.f; see Genz (2004) for the background methodology. A similar fact takes place when df=2 with function biv.nt.prob; note however that the underlying Fortran code is taken from mvtdstpack.f, not from tvpack.f. If pmt is called with d>3, this is converted into a suitable call to sadmvt.

If sqrt=NULL (default value), the working of rmt involves computation of a square root of S via the Cholesky decomposition. If a non-NULL value of sqrt is supplied, it is assumed that it represents a square root of the scale matrix, otherwise represented by S, whose value is ignored in this case. This mechanism is intended primarily for use in a sequence of calls to rmt, all sampling from a distribution with fixed scale matrix; a suitable matrix sqrt can then be computed only once beforehand, avoiding that the same operation is repeated multiple times along the sequence of calls. For examples of use of this argument, see those in the documentation of rmnorm. Another use of sqrt is to supply a different form of square root of the scale matrix, in place of the Cholesky factor.

For efficiency reasons, rmt does not perform checks on the supplied arguments.

Value

dmt returns a vector of density values (possibly log-transformed); pmt and sadmvt return a single probability with attributes giving details on the achieved accuracy, provided x of pmnorm is a vector; rmt returns a matrix of n rows of random vectors, or a vector in case n=1 or d=1.

References

Genz, A.: Fortran 77 code in files mvt.f, mvtdstpack.f and tvpack, downloaded in 2005 and again in 2007 from his software webpage, whose URL as of 2020-06-01 was https://www.math.wsu.edu/faculty/genz/software/software.html

Genz, A. (2004). Numerical computation of rectangular bivariate and trivariate normal and t probabilities. Statistics and Computing 14, 251-260.

Dunnett, C.W. and Sobel, M. (1954). A bivariate generalization of Student's t-distribution with tables for certain special cases. Biometrika 41, 153–169.

Author

FORTRAN 77 code of SADMVT, MVTDSTPACK, TVPACK and many auxiliary functions by Alan Genz; some additional auxiliary functions by people referred to within his programs; interface to R and additional R code (for dmt, rmt etc.) by Adelchi Azzalini.

Note

The attributes error and status of the probability returned by sadmvt and by pmt (the latter only if x is a vector and d>2) indicate whether the function had a normal termination, achieving the required accuracy. If this is not the case, re-run the function with a higher value of maxpts.

Examples

x <- seq(-2,4,length=21)
y <- 2*x+10
z <- x+cos(y) 
mu <- c(1,12,2)
Sigma <- matrix(c(1,2,0,2,5,0.5,0,0.5,3), 3, 3)
df <- 4
f  <- dmt(cbind(x,y,z), mu, Sigma,df)
p1 <- pmt(c(2,11,3), mu, Sigma, df)
p2 <- pmt(c(2,11,3), mu, Sigma, df, maxpts=10000, abseps=1e-8)
x  <- rmt(10, mu, Sigma, df)
p  <- sadmvt(df, lower=c(2,11,3), upper=rep(Inf,3), mu, Sigma) # upper tail
#
p0 <- pmt(c(2,11), mu[1:2], Sigma[1:2,1:2], df=5)
p1 <- biv.nt.prob(5, lower=rep(-Inf,2), upper=c(2, 11), mu[1:2], Sigma[1:2,1:2])
p2 <- sadmvt(5, lower=rep(-Inf,2), upper=c(2, 11), mu[1:2], Sigma[1:2,1:2]) 
c(p0, p1, p2, p0-p1, p0-p2)
#> [1] 0.3356944 0.3356944 0.3356944 0.0000000 0.0000000