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, 2001).

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

_V_a_l_u_e:

     The distribution function is returned with attributes 

   error: estimated absolute error and

     msg: status messages.

_A_u_t_h_o_r(_s):

     Fortran Code by Alan Genz <AlanGenz@wsu.edu> and Frank Bretz
     <bretz@ifgb.uni-hannover.de>, R port by Torsten Hothorn
     <Torsten.Hothorn@rzmail.uni-erlangen.de>

_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. (2001), Methods for the computation of
     multivariate t-probabilities. (submitted)

     Edwards D. and Berry, Jack J. (1987), The efficiency of
     simulation-based  multiple comparisons. Biometrics, 43, 913-928

_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)
     cv <- C %*% V %*% t(C)
     cr <- matrix(rep(0, ncol(cv)^2), ncol=ncol(cv))
     for (i in 1:5) {
       for (j in 1:5) {
         cr[i,j] <- cv[i,j]/sqrt(cv[i,i]*cv[j,j] )
       }
     }
     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)))

