fastSP: Functions for fast calculation of stratification pattern...

In bertcarnell/SOAs: Creation of Stratum Orthogonal Arrays

Functions for fast calculation of stratification pattern according to Tian and Xu 2023

Description

Functions for fast calculation of stratification pattern according to Tian and Xu 2023

Usage

fastSP(D, s, maxwt = NULL, K = NULL, y0 = NULL, tol = 1e-05)

Arguments

D	design with number of levels a power of s
s	prime or prime power on which D is based
maxwt	integer number; maximum weight for which the pattern is to be calculated
K	integer number of summands. Can also be "max" for indicating that all summands are requested (Theorem 3 of Tian and Xu 2023+). In Theorem 4 of Tian and Xu (2023+), larger K provide better accuracy. If maxwt=NULL, the default (NULL) is that all summands are requested (i.e., Theorem 3 of Tian and Xu 2023+). Otherwise, the default is the maximum of maxwt and a default based on y0 according to a formula by Tian and Xu (2023+) (yields smaller K for smaller y0).
y0	small number that drives accuracy. The default (NULL) uses 1/s for maximum K and y0=0.1 for smaller values of K. See the Details section for a discussion of this parameter.
tol	tolerance for checking whether the imaginary part is zero

Details

The function was modified from the code provided with Tian and Xu (2023+).

Per default (maxwt=NULL and K=NULL), or when the user chooses K="max" in spite of specifying a value for maxwt, fastSP calculates the entire stratification pattern and uses all summands of Tian and Xu's (2023+) Equation (4) of their Theorem 3,
with one important modification: y_j is not chosen as the j-th power of the m*el-th complex rooth of the unity but as that complex number divided by s, as this appears to yield numerically stabler results. (It neutralizes the large powers of s that arise by multiplying the weighted similarities of Lemma 1 for obtaining the summands in E(D,y) in formula (3)). Limited experimentation with different values of s showed that this approach did indeed yield reasonably stable results, for example for the SOA(64,20,8,3) for which the version with unmodified powers of roots of the unity ran into problems.
It is straightforward to verify that Theorem 3 remains valid for this modified choice of y_j.

It is possible and often advisable to calculate only a smaller number of entries, for saving resources and also because the later entries are less interesting and less accurate. If the argument maxwt is specified and K is left unspecified (K=NULL), K is calculated as the maximum of Tian and Xu's (2023+) proposed default for their approximation formula (6), depending on y0, which is set to 0.1, if also unspecified. Tian and Xu (2023+) recommended to use y0 values between 0.001 and 0.1, when using formula (6) of Theorem 4.

For obtaining the original behavior of Tian and Xu's (2023+) implementation of their Theorem 3 (not desirable for larger situations), choose y0=1. Note that y0=1 with a specified maxwt and K=NULL yields an error, because the default formula for K does not work for y0=1.

Also consider the Note section regarding numerical considerations.

Value

an object of class fastSP, with attributes call, K, y0, and possibly message. The object itself is a stratification pattern or the first maxwt elements of the stratification pattern (default: all elements).
If K is less than the maximum length of the stratification pattern (Kmax element of attribute K) the returned values are approximations (more accurate for larger K). If the object has a message attribute, this attribute indicates which positions of the pattern must be considered as problematic because the imaginary part was non-zero.

Note

Even the exact pattern (obtained with maximum K) must be considered with caution because of potential numerical problems. Often, the creation process of a GSOA implies that the first few elements are zeroes. If this is the case, the degree of inaccuracy may be assessed from these elements. Furthermore, warnings of non-zero imaginary parts indicate similar problems. If unsure about the accuracy, it may also be an option to use function Spattern with a small maxwt argument (for resource reasons) in order to obtain exact values for the first very few entries of the stratification pattern.aus

References

For full detail, see SOAs-package.

Tian, Y. and Xu, H. (2023+)

Examples

## SOA(32,9,8,3) from Shi and Tang (2020)
soa32x9 <- t(matrix(c(7,3,6,2,7,3,6,2,4,0,5,1,4,0,5,1,5,1,4,0,5,1,4,0,6,2,7,3,6,2,7,3,
                  7,7,2,2,5,5,0,0,6,6,3,3,4,4,1,1,5,5,0,0,7,7,2,2,4,4,1,1,6,6,3,3,
                  7,5,6,4,3,1,2,0,4,6,5,7,0,2,1,3,7,5,6,4,3,1,2,0,4,6,5,7,0,2,1,3,
                  7,7,4,4,5,5,6,6,2,2,1,1,0,0,3,3,7,7,4,4,5,5,6,6,2,2,1,1,0,0,3,3,
                  7,5,6,4,5,7,4,6,6,4,7,5,4,6,5,7,3,1,2,0,1,3,0,2,2,0,3,1,0,2,1,3,
                  7,1,0,6,3,5,4,2,4,2,3,5,0,6,7,1,5,3,2,4,1,7,6,0,6,0,1,7,2,4,5,3,
                  7,1,2,4,7,1,2,4,2,4,7,1,2,4,7,1,5,3,0,6,5,3,0,6,0,6,5,3,0,6,5,3,
                  7,3,2,6,5,1,0,4,4,0,1,5,6,2,3,7,3,7,6,2,1,5,4,0,0,4,5,1,2,6,7,3,
                  7,1,4,2,3,5,0,6,2,4,1,7,6,0,5,3,3,5,0,6,7,1,4,2,6,0,5,3,2,4,1,7
), nrow=9, byrow = TRUE))

## complete pattern according to theorem 3
## (y0=1/2 or y0=1 does not make a difference for this small example)
a <- fastSP(soa32x9, 2); round(a,7)
a <- fastSP(soa32x9, 2, y0=1); round(a,7)

## only the first five positions (K=9 calculated and used based on y0=0.1)
## not very accurate
a <- fastSP(soa32x9, 2, maxwt=5); round(a,7)
## more accurate (K=5 used, based on y0=0.01)
a <- fastSP(soa32x9, 2, maxwt=5, y=0.01); round(a,7)
## even more accurate (K=9 used with y0=0.01)
a <- fastSP(soa32x9, 2, maxwt=5, y=0.01, K=9); round(a,7)

# example code

bertcarnell/SOAs documentation built on Nov. 21, 2023, 11:25 a.m.