pmvt                 package:mvtnorm                 R Documentation

_M_u_l_t_i_v_a_r_i_a_t_e _t _D_i_s_t_r_i_b_u_t_i_o_n

_D_e_s_c_r_i_p_t_i_o_n:

     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.

_U_s_a_g_e:

     pmvt(lower=-Inf, upper=Inf, delta=rep(0, length(lower)),
          df=1, corr=NULL, sigma=NULL, maxpts = 25000, abseps = 0.001,
          releps = 0)

_A_r_g_u_m_e_n_t_s:

   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. 

_D_e_t_a_i_l_s:

     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)(2pi)^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  <URL:
     http://www.bioinf.uni-hannover.de/~bretz/>  and <URL:
     http://www.sci.wsu.edu/math/faculty/genz/homepage>.

     If 'df = 0', normal probabilities are returned.

_V_a_l_u_e:

     The evaluated distribution function is returned with attributes 

   error: estimated absolute error and

     msg: status messages.

_R_e_f_e_r_e_n_c_e_s:

     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.

_S_e_e _A_l_s_o:

     'qmvt'

_E_x_a_m_p_l_e_s:

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