The probability density function, the distribution function and random number 
  generation for a d-dimensional Student's t random variable.
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)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.
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’.
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.
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).
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.
a logical value(default value is FALSE); if TRUE, 
     the logarithm of the density is computed.
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.
the number of random vectors to be generated
a numeric vector of lower integration limits of 
     the density function; must be of maximal length 20; 
     +Inf and -Inf entries are allowed.
a numeric vector of upper integration limits 
     of the density function; must be of maximal length 20; 
     +Inf and -Inf entries are allowed
the maximum number of function evaluations 
               (default value: 2000*d)
absolute error tolerance (default value: 1e-6).
relative error tolerance (default value: 0).
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.
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, available 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.
Genz, A.:  Fortran 77 code in files mvt.f, mvtdstpack.f  
  and codetvpack, downloaded in 2005 and again in 2007 from his webpage, 
  whose URL as of 2020-06-01 is 
  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.
dt, 
       rmnorm for use of argument sqrt,
       plot_fxy for plotting 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)
Run the code above in your browser using DataLab