rankSTB {STB}

Rank-Based Algorithm for Computing 100(1-alpha)% Simulataneous Tolerance Bounds.

Description

Implementation of a rank-based algorithm for constructing 100(1-alpha)% STBs as outlined in the reference.

Usage

rankSTB(mat, alpha = 0.05)

Arguments

mat

(numeric) matrix of dimension (N, n), where i-th row corrsponds to ordered values of the i-th simulation

alpha

(numeric) value defining the desired coverage as 100(1-alpha)%

Details

This function is a performance optimized version of the original rank-based algorithm avoiding the time-consuming iteration. In principle it sorts out simulation results which have at least one extreme order statistic untill exactly 100(1-alpha)% of all simulation results remain. From these, bounds of the STB are constructed determining extreme-values per order-statistic (column).

This implementation also corrects step 4) of the published algorithm, which has to find those indices of elements being equal to "min(c)" OR being equal to "N-min(c)+1". This reflects the construction of vector "c", where max. rank-values are transformed to min. rank-values. In step 6) the "N_k-1-(1-alpha)*N largest" elements of "d_l^theta" have to be selected, which needs also correction.

Parallel processing did not minimize the computation time in contrast to the algorithms for computing the quantile-based algorithm. Thus, parallel processing is not supported for this algorithm.

Value

(list) with two elements:

Q

(matrix) 1st row stores lower bounds, 2nd row upper bounds

cov

(numeric) value corresponding the coverage

Author(s)

Andre Schuetzenmeister [e-mail]

References

Schuetzenmeister, A. and Piepho, H.P. (2012). Residual analysis of linear mixed models using a simulation approach. Computational Statistics and Data Analysis, 56, 1405-1416

Examples

## Not run: 
# for following problem size the rank-based algo
# outperforms the quantile based one, although,
# ran serially
mat <- matrix(rnorm(10000*100), ncol=100)
mat <- t(apply(mat, 1, sort))
system.time(stb.rank <- rankSTB(mat))
system.time(stb.q.R  <- getSTB(mat))
system.time(stb.q.C  <- fastSTB(mat))
x <- apply(mat, 2, mean)
plot(x,x, ylim=c(-5,5))
lines(x, stb.q.R$Q[1,], col="blue", lwd=2)
lines(x, stb.q.R$Q[2,], col="blue", lwd=2)
lines(x, stb.q.C$Q[1,], col="red",  lwd=2)
lines(x, stb.q.C$Q[2,], col="red",  lwd=2)
lines(x, stb.rank$Q[1,],  col="cyan", lwd=2)
lines(x, stb.rank$Q[2,],  col="cyan", lwd=2)
legend("top", legend=c("R-quantile", "C-quantile", "rank-based"), 
       fill=c("blue", "red", "cyan"))

# varying Ncpu for the C-implementation of the quantile-based algo
system.time(stb.q.C  <- fastSTB(mat, Ncpu=4))
system.time(stb.q.C  <- fastSTB(mat, Ncpu=6))
system.time(stb.q.C  <- fastSTB(mat, Ncpu=8))
system.time(stb.q.C  <- fastSTB(mat, Ncpu=10))
system.time(stb.q.C  <- fastSTB(mat, Ncpu=12))

## End(Not run)