SOAs (version 1.4-1)

fastSP: Function for fast calculation of stratification pattern according to Tian and Xu 2023

Description

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

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.

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 \((S_1(D),...,S_{m\ell}(D))\) of the design \(D\) with \(m\) columns in \(s^\ell\) levels based on solving the system of \(m\ell\) equations $$E(D;y_i) - 1 = \sum_{j=1}^{m\ell}y_i^{j}S_j(D)$$ with \(E(D;y_i)\) as defined in Tian and Xu's (2023+) equation (3) and \(y_i=1/s\cdot \tilde\omega_i\) with \(\tilde\omega_i\) the \(m\ell\)th complex root of the unity. This system arises from Tian and Xu's (2023+) Theorem 1, and its solution is stated in their equation (4) in Theorem 3 for a numerically less fortunate choice of \(y_i\) values. (Theorem 3 of Tian and Xu 2023+ uses \(\tilde\omega_i\) instead of the above \(y_i\) values; \(y_i=\tilde\omega_i/s\) improves the numerical stability by neutralizing large powers of s that otherwise arise in the summands of \(E(D;y_i)\); Example 6 of the paper that had numerical problems for Theorem 3 worked well for the entire pattern with this implementation.) The implementation for K="all" can also be considered a special case of Tian and Xu's (2023+) Theorem 4 with \(z=1/s\) and \(K=m\ell\). (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.)

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.

Consider also the Note section regarding recommendations for accuracy checks.

References

For full detail, see SOAs-package.

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

See Also

Spattern() for exact calculations of the first few elements of the stratification pattern for small to moderate situations. That function additionally provides a dimension by weight table.

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
a <- fastSP(soa32x9, 2); round(a,7)

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