pmvss {mvgb} | R Documentation |
Multivariate Subgaussian Stable Distribution Probabilities (via positive stable variates)
Description
Computes the the distribution function of the multivariate subgaussian stable
distribution for arbitrary limits, alpha,
shape matrices, and location vectors.
This function is unlike mvtnorm::pmvt
in two ways:
The QRVSN method is used for every dimension
n
(including bivariate and trivariate),The QRVSN method on positive stable variates, not chi/sqrt(nu).
Usage
pmvss(
lower,
upper,
alpha,
Q,
delta = rep(0, NROW(Q)),
maxpts = 25000,
abseps = 0.001,
releps = 0
)
Arguments
lower |
lower bounds of integration must be length n, finite |
upper |
upper bounds of integration must be length n, finite |
alpha |
real number between 0 and 2 that cannot have more than 6 digits precision. 0.123456 is okay; 0.1234567 is not. Thus alpha in [0.000001, 1.999999]. |
Q |
shape matrix |
delta |
location vector must have length equal to the number of rows of Q. Defaults to the 0 vector. |
maxpts |
(description from FORTRAN code) INTEGER, maximum number of function values allowed. This parameter can be used to limit the time. A sensible strategy is to start with MAXPTS = 1000*N, and then increase MAXPTS if ERROR is too large. (description from mvtnorm::GenzBretz) maximum number of function values as integer. The internal FORTRAN code always uses a minimum number depending on the dimension. (for example 752 for three-dimensional problems). |
abseps |
absolute error tolerance |
releps |
relative error tolerance as double. |
Value
Returns the variables from the MVTDST function (QRVSN algorithm):
N INTEGER, the number of variables. NU INTEGER, the number of degrees of freedom. If NU < 1, then an MVN probability is computed. LOWER DOUBLE PRECISION, array of lower integration limits. UPPER DOUBLE PRECISION, array of upper integration limits. INFIN INTEGER, array of integration limits flags: if INFIN(I) < 0, Ith limits are (-infinity, infinity); if INFIN(I) = 0, Ith limits are (-infinity, UPPER(I)]; if INFIN(I) = 1, Ith limits are [LOWER(I), infinity); if INFIN(I) = 2, Ith limits are [LOWER(I), UPPER(I)]. CORREL DOUBLE PRECISION, array of correlation coefficients; the correlation coefficient in row I column J of the correlation matrix should be stored in CORREL( J + ((I-2)*(I-1))/2 ), for J < I. The correlation matrix must be positive semi-definite. DELTA DOUBLE PRECISION, array of non-centrality parameters. MAXPTS INTEGER, maximum number of function values allowed. This parameter can be used to limit the time. A sensible strategy is to start with MAXPTS = 1000*N, and then increase MAXPTS if ERROR is too large. ABSEPS DOUBLE PRECISION absolute error tolerance. RELEPS DOUBLE PRECISION relative error tolerance. error DOUBLE PRECISION estimated absolute error, with 99% confidence level. value DOUBLE PRECISION estimated value for the integral inform INTEGER, termination status parameter: if INFORM = 0, normal completion with ERROR < EPS; if INFORM = 1, completion with ERROR > EPS and MAXPTS function values used; increase MAXPTS to decrease ERROR; if INFORM = 2, N > 1000 or N < 1. if INFORM = 3, correlation matrix not positive semi-definite. RND ignore; this initializes RNG
References
Genz, A. and Bretz, F. (2002), Methods for the computation of multivariate t-probabilities. Journal of Computational and Graphical Statistics, 11, 950–971.
http://www.math.wsu.edu/faculty/genz/homepage
http://www.math.wsu.edu/faculty/genz/software/fort77/mvtdstpack.f
Examples
Q = structure(c(1, 0.85, 0.85, 0.85, 0.85,
0.85, 1 , 0.85, 0.85, 0.85,
0.85, 0.85 , 1 , 0.85, 0.85,
0.85, 0.85 , 0.85 , 1 , 0.85,
0.85, 0.85 , 0.85 , 0.85 , 1 ),
.Dim = c(5L,5L))
## default maxpts=25000 doesn't finish with error < abseps
mvgb::pmvss(lower=rep(-1,5),
upper=rep(1,5),
alpha=1,
Q=Q,
maxpts=25000)[c("value","inform","error")]
## increase maxpts to get inform value 0 (that is, error < abseps)
mvgb::pmvss(lower=rep(-1,5),
upper=rep(1,5),
alpha=1,
Q=Q,
maxpts=25000*350)[c("value","inform","error")]
set.seed(10)
shape_matrix <- structure(c(1, 0.9, 0.9, 0.9, 0.9, 0.9, 1, 0.9, 0.9, 0.9, 0.9,
0.9, 1, 0.9, 0.9, 0.9, 0.9, 0.9, 1, 0.9, 0.9, 0.9, 0.9, 0.9,
1), .Dim = c(5L, 5L))
mvgb::pmvss(lower=rep(-2,5),
upper=rep(2,5),
alpha=1.7,
Q=shape_matrix,
delta=rep(0,5),
maxpts=25000*30)[c("value","inform","error")]