Package {HDElliptical}

HDElliptical: high-dimensional elliptical methods

Description

HDElliptical provides R and Rcpp implementations accompanying High-Dimensional Data Analysis for Elliptically Symmetric Distributions. Functions use observations in rows and variables in columns throughout.

Author(s)

Maintainer: Long Feng flnankai@nankai.edu.cn (Copyright holder for the HDElliptical implementation) [copyright holder]

Other contributors:

• Dan Zhuang (SEMC methodology and the external GEMcluster reference implementation used for validation; no GEMcluster source code is included) [contributor]

See Also

Useful links:

• https://github.com/flnankai/HDElliptical

• Report bugs at https://github.com/flnankai/HDElliptical/issues

Angular central Gaussian log-likelihood

Description

Evaluates the angular central Gaussian log-likelihood for nonzero row vectors. Row lengths need not equal one because they cancel from the quadratic contribution.

Usage

acg_loglik(x, shape, include_constant = FALSE)

Arguments

Value

A scalar log-likelihood.

References

Tyler, D. E. (1987). A distribution-free M-estimator of multivariate scatter. Annals of Statistics, 15, 234–251. doi:10.1214/aos/1176350263.

Examples

directions <- rbind(c(1, 0), c(0, 1), c(1, 1))
acg_loglik(directions, diag(2))

Original Bai–Saranadasa two-sample mean test

Description

Tests equality of two multivariate means under a common covariance matrix using the original Bai–Saranadasa statistic. This is deliberately a separate API from covariance-unrestricted Chen–Qin tests: its numerator is

\xi = \|\bar X_1 - \bar X_2\|^2 - \frac{N}{n_1 n_2}\operatorname{tr}(S_p),

where S_p is the unbiased pooled covariance matrix and N=n_1+n_2.

Usage

bai_saranadasa_two_sample_test(x, y)

Arguments

Details

The traces of S_p and S_p^2 are computed using the smaller of a primal p \times p matrix and a dual N \times N Gram matrix. Both groups are divided by the same global numerical scale before any moments are formed. The standardised statistic is computed on that safe scale, while reported means, traces, numerator, and variance are converted back according to their physical dimensions. Consequently a reported fourth-order component can be infinite when its value exceeds double range, or zero when it is below the representable range, without making the standardised statistic or p-value non-finite. The finite internal normalising variance is retained in diagnostics. The estimated null variance must be strictly positive; no ridge, absolute-value repair, or variance flooring is applied.

Value

An object of class c("hd_location_test", "htest"). Its components include the original numerator xi, pooled covariance traces, the bias-corrected estimate of \operatorname{tr}(\Sigma^2), and the estimated null variance.

References

Bai, Z. and Saranadasa, H. (1996). Effect of high dimension: by an example of a two sample problem. Statistica Sinica, 6, 311–329.

Examples

set.seed(22)
x <- matrix(rnorm(24), nrow = 8, ncol = 3)
y <- matrix(rnorm(27, mean = 0.2), nrow = 9, ncol = 3)
bai_saranadasa_two_sample_test(x, y)

BASIC bias-adjusted spatial-sign shape estimator

Description

Applies the real-valued BASIC approximate inverse map independently to the eigenvalues of the fitted SSCM. For a shape eigenvalue \lambda, the map used by Raninen and Ollila is

\widetilde\delta(\lambda)=\frac12\int_0^1 \frac{\lambda}{1-t+t\lambda}(1-t)^{p/2-1}\,dt.

The implementation evaluates the equivalent smooth integral after 1-t=v^2 with paired Gauss–Legendre rules, brackets each inverse root, and uses bisection. It is independent of the official lookup table: no interpolation, spline extrapolation, eigenvalue floor, or out-of-range repair is used. Raw inverted eigenvalues are finally normalized to sum p.

Usage

basic_shape(
  x,
  center = c("spatial", "mean", "none"),
  tol = 1e-08,
  max_iter = 1000L,
  strict = TRUE,
  keep_signs = FALSE,
  quadrature_order = 64L,
  inversion_tol = 1e-10,
  integration_tol = 1e-10,
  inversion_max_iter = 100L,
  max_bracket_iter = 100L
)

Arguments

Details

Boundary SSCM eigenvalues equal to zero map exactly to zero and are reported. For p = 1, shape is identically one and no numerical inversion is needed.

Value

An hd_shape_estimator with the BASIC shape, SSCM eigensystem, raw and normalized inverse eigenvalues, mapped values, integration and inversion errors, brackets, iterations, and boundary diagnostics.

References

Raninen, E. and Ollila, E. (2022). Bias adjusted sign covariance matrix. IEEE Signal Processing Letters, 29, 339–343. doi:10.1109/LSP.2021.3134940.

Examples

x <- matrix(c(-2, 0, -1, 2, 0, -2, 1, 1, 2, -1, 3, 2),
            ncol = 2, byrow = TRUE)
basic_shape(x, quadrature_order = 32)

BASICS shrinkage plus bias-adjusted shape estimator

Description

Computes the RSSCM data weight, forms \widehat\Lambda_{RSSCM}=\widehat\alpha pS_{sgn}+ (1-\widehat\alpha)I_p, and applies the same self-contained BASIC inverse map as basic_shape() to the eigenvalues of \widehat\Lambda_{RSSCM}/p. This follows the BASICS construction while avoiding the official code's spline extrapolation. The RSSCM and BASIC numerical diagnostics are both retained.

Usage

basics_shape(
  x,
  center = c("spatial", "mean", "none"),
  tol = 1e-08,
  max_iter = 1000L,
  strict = TRUE,
  keep_signs = FALSE,
  quadrature_order = 64L,
  inversion_tol = 1e-10,
  integration_tol = 1e-10,
  inversion_max_iter = 100L,
  max_bracket_iter = 100L
)

Arguments

Value

An hd_shape_estimator containing the BASICS shape, RSSCM, its raw/projected shrinkage weight, inverse-map details, and full diagnostics.

References

Raninen, E. and Ollila, E. (2022). Bias adjusted sign covariance matrix. IEEE Signal Processing Letters, 29, 339–343. doi:10.1109/LSP.2021.3134940.

Examples

x <- matrix(c(-2, 0, -1, 2, 0, -2, 1, 1, 2, -1, 3, 2),
            ncol = 2, byrow = TRUE)
basics_shape(x, quadrature_order = 32)

Bickel–Levina hard-thresholded covariance estimator

Description

Forms the centered sample covariance with divisor n, as in Bickel and Levina (2008), and applies

s_{ij}\,1\{|s_{ij}|>\lambda\}.

When threshold = NULL, \lambda=C\sqrt{\log(p)/n}. The book draft first defines a divisor-n-1 covariance and then cites the primary divisor-n formula. divisor makes this finite-sample distinction explicit; the default is the primary-paper convention. Thresholding need not preserve positive definiteness and no eigenvalue repair is performed.

Usage

bickel_levina_covariance_threshold(
  x,
  threshold = NULL,
  constant = 1,
  center = TRUE,
  divisor = c("n", "n-1"),
  threshold_diagonal = TRUE
)

Arguments

Value

An hd_covariance_estimator list whose estimate field is the thresholded covariance matrix.

References

Bickel, P. J. and Levina, E. (2008). Annals of Statistics, 36, 2577–2604.

Examples

set.seed(31)
x <- matrix(rnorm(40), nrow = 10, ncol = 4)
bickel_levina_covariance_threshold(x, threshold = 0.2)

Book Gaussian alpha Cauchy benchmark

Description

Implements the formula-complete Gaussian benchmark printed in Chapter 4, kept distinct from gaussian_alpha_combination_test(), which implements the primary Feng–Lan–Liu–Ma Bonferroni rule. With unrestricted OLS residual variances \hat\sigma_i^2, residual correlation \hat R, and intercept estimates \hat\alpha, this function uses

T_{sum} = \frac{T\sum_i\hat\alpha_i^2/\hat\sigma_i^2-N} {\{2\mathrm{tr}(\hat R^2)\}^{1/2}}

and the maximum squared OLS intercept t statistic. Their upper-tail normal and extreme-value p-values are combined by the equal-weight Cauchy rule.

Usage

book_gaussian_alpha_cauchy_test(returns, trace_R2, factors = NULL)

Arguments

Details

The manuscript denotes the trace term abstractly and does not prescribe a unique finite-sample estimator for this Gaussian benchmark. Consequently, trace_R2 is required rather than silently replacing it by a residual plug-in or a bias-corrected alternative.

Value

An upper-tail alpha-test object explicitly labelled as a book benchmark rather than a primary named procedure.

References

High-Dimensional Data Analysis for Elliptically Symmetric Distributions, Chapter 4, equations for the Gaussian alpha benchmark.

Examples

f <- cbind(seq(-1, 1, length.out = 12))
y <- outer(seq_len(12), 1:3, function(i, j) cos(i / 2 + j))
book_gaussian_alpha_cauchy_test(y, trace_R2 = 3, factors = f)

Cai–Liu adaptive covariance thresholding

Description

The primary-paper pilot covariance and variability estimate are

\hat\sigma_{ij}=n^{-1}\sum_k Z_{ki}Z_{kj},\qquad \hat\theta_{ij}=n^{-1}\sum_k(Z_{ki}Z_{kj}-\hat\sigma_{ij})^2,

with entrywise threshold

\lambda_{ij}=\delta \sqrt{\hat\theta_{ij}\log(p)/n}.

Both quantities therefore use the same divisor-n covariance by default. Selecting divisor = "n-1" is explicit and is reported in the diagnostics; it is not the primary finite-sample formula.

Usage

cai_liu_adaptive_covariance_threshold(
  x,
  delta = 2,
  rule = c("hard", "soft", "scad", "adaptive_lasso"),
  scad_a = 3.7,
  adaptive_eta = 1,
  center = TRUE,
  divisor = c("n", "n-1"),
  threshold_diagonal = TRUE
)

Arguments

Value

An hd_covariance_estimator object retaining \hat\theta and every entrywise threshold.

References

Cai, T. and Liu, W. (2011). Journal of the American Statistical Association, 106, 672–684. arXiv:1102.2237.

Examples

set.seed(33)
x <- matrix(rnorm(40), nrow = 10, ncol = 4)
cai_liu_adaptive_covariance_threshold(x, delta = 2)

Cai–Liu–Xia precision-adjusted two-sample maximum test

Description

Tests equality of two high-dimensional mean vectors using the precision-adjusted maximum statistic of Cai, Liu, and Xia (2014). The calibration is the type-I extreme-value limit

F(g)=\exp\{-\pi^{-1/2}\exp(-g/2)\},

for

G=M-2\log(p)+\log\{\log(p)\}.

Usage

cai_liu_xia_two_sample_test(
  x,
  y,
  precision = NULL,
  precision_source = c("adaptive", "oracle", "feasible"),
  threshold_delta = 2,
  eigen_floor = sqrt(.Machine$double.eps)
)

Arguments

Details

precision_source deliberately separates three statistically different paths. With "oracle", precision is a known population precision matrix and the denominator is its diagonal. With "feasible", precision is a user-supplied estimate and the denominator is formed from empirical within-group variances of the transformed observations, each with divisor n_k, as in equations (6)–(7) of the paper. With "adaptive", the function estimates a precision matrix by hard-thresholding the pooled covariance and then uses that same feasible denominator. It never uses diag(precision) to standardise an estimated-precision statistic.

For the adaptive backend, the pooled covariance and \widehat\theta_{ij} both use divisor N=n_1+n_2, and

\lambda_{ij}=\delta \sqrt{\widehat\theta_{ij}\log(p)/N}.

If the thresholded covariance is not numerically positive definite, its eigenvalues are raised to eigen_floor times a data-scale reference before inversion. The floor, the unadjusted and adjusted extreme eigenvalues, and the number of modified eigenvalues are all returned in diagnostics. Both samples are centered relative to observed anchors before precision transformation, avoiding cancellation under a large common translation. The adaptive covariance and fourth-order threshold moments are also formed on an internal common scale and then mapped back to the original units.

Value

An object of class c("hd_location_test", "htest"). The reported statistic G has the extreme-value calibration, while raw.statistic contains M. Coordinatewise transformed scores, their denominators, the maximizing coordinate index and name, and the precision matrix used by the test are retained in components. diagnostics identifies the precision path, denominator construction, and every adaptive-threshold eigenvalue adjustment.

References

Cai, T. T., Liu, W., and Xia, Y. (2014). Two-sample test of high dimensional means under dependence. Journal of the Royal Statistical Society: Series B, 76, 349–372.

Cai, T. T. and Liu, W. (2011). Adaptive thresholding for sparse covariance matrix estimation. Journal of the American Statistical Association, 106, 672–684.

Examples

set.seed(29)
x <- matrix(rnorm(80), 20, 4)
y <- matrix(rnorm(96, 0.2), 24, 4)
cai_liu_xia_two_sample_test(x, y)

omega <- diag(4)
cai_liu_xia_two_sample_test(
  x, y, precision = omega, precision_source = "oracle"
)

Bartlett likelihood-ratio approximation for canonical correlations

Description

Tests the tail null that canonical correlations after null_rank are zero. For null_rank = k, the statistic uses

\Lambda_k=\prod_{j=k+1}^{m}(1-\widehat\rho_j^2),\qquad -\{n-1-(p_x+p_y+1)/2\}\log\Lambda_k,

with chi-squared degrees of freedom (p_x-k)(p_y-k). The special case k = 0 is the full independence test printed in the book.

Usage

cca_bartlett_test(object, null_rank = 0L)

Arguments

Value

An object of class htest with Wilks' Lambda and the Bartlett approximation diagnostics.

References

Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis, 3rd ed. Wiley.

Examples

basis <- stats::poly(seq_len(12), degree = 4)
x <- basis[, 1:2, drop = FALSE]
y <- cbind(0.8 * basis[, 1] + 0.6 * basis[, 3],
           0.3 * basis[, 2] + sqrt(0.91) * basis[, 4])
colnames(x) <- c("x1", "x2")
colnames(y) <- c("y1", "y2")
fit <- classical_cca(x, y)
cca_bartlett_test(fit, null_rank = 0)

Coordinate CUSUMs and Rice difference variances

Description

Computes the literal gamma=0 and gamma=1/2 CUSUM arrays together with sum(diff(x)^2)/(2(n-1)). Columnwise anchors and long-double accumulators reduce cancellation under a common translation.

Usage

ch4_cp_cusum_cpp(x)

Arguments

Value

A list of two CUSUM matrices and the difference variance.

Wang–Feng finite-difference moment estimators

Description

Implements the displayed leave-four and leave-three estimators in the primary paper. The exclusion set is literally ⁠{2,...,n} \\ {i1,...,im}⁠ and no bridging difference is inserted.

Usage

ch4_cp_dms_moments_cpp(x)

Arguments

Value

Trace, radial fourth-moment, and term-level diagnostics.

ERHT companion center and variance

Description

Computes the ordered off-diagonal expression directly, avoiding subtraction of two nearly equal non-negative matrix products.

Usage

ch4_cp_erht_moments_cpp(companion, beta)

Arguments

Value

Kappa, sigma squared, and the ordered off-diagonal sum.

Ordered-pair squared inner-product sum

Description

Ordered-pair squared inner-product sum

Usage

ch4_cp_ordered_pair_square_sum_cpp(signs)

Arguments

Value

⁠sum_{i != j} (u_i' u_j)^2⁠.

Joint scaled spatial-median and diagonal fixed point

Description

Implements the three updates stated by Liu, Feng, Peng, and Wang. The common diagonal scale inherited from sample-variance initialization is retained; only relative diagonal scale is identified.

Usage

ch4_cp_scaled_hr_cpp(x, tol = 1e-08, max_iter = 1000L, zero_tol = 0)

Arguments

Value

Fixed point and convergence diagnostics.

Spatial median with a coincident-point certificate

Description

Uses the modified Weiszfeld map. At an iterate coinciding with observations, the exact subgradient condition is checked rather than adding jitter.

Usage

ch4_cp_spatial_median_cpp(x, tol = 1e-08, max_iter = 1000L, zero_tol = 0)

Arguments

Value

Median and convergence diagnostics.

Chen–Qin two-sample high-dimensional mean test

Description

Tests equality of two mean vectors with the covariance-unrestricted U-statistic of Chen and Qin (2010). Its numerator is

T_{CQ}=\|\bar X-\bar Y\|^2- \operatorname{tr}(S_X)/n_1-\operatorname{tr}(S_Y)/n_2.

The variance estimator is the original paper's leave-out estimator of the two within-sample covariance traces and their cross trace. This is distinct from the common-covariance Bai–Saranadasa test.

Usage

chen_qin_two_sample_test(x, y)

Arguments

Details

The numerator is exactly invariant to a common translation. The original finite-sample leave-out variance estimator, however, is not translation invariant; this known property is retained rather than silently replacing the published estimator. Both samples are internally divided by one common global scale. The returned components include both the original-unit quantities and their finite internally scaled counterparts. An original-unit fourth-order quantity can legitimately overflow to Inf or underflow to zero when the input itself spans the limits of double precision; the reported statistic is reconstructed from the scaled quantities and remains well defined. No ridge, absolute-value repair, or variance flooring is applied.

Value

An object of class c("hd_location_test", "htest"). The components field includes the U-statistic numerator, leave-out trace estimates A1, A2, and A12, the estimated null variance, and their .scaled counterparts used for stable calibration.

References

Chen, S. X. and Qin, Y.-L. (2010). A two-sample test for high-dimensional data with applications to gene-set testing. Annals of Statistics, 38, 808–835.

Examples

set.seed(24)
x <- matrix(rnorm(21), nrow = 7, ncol = 3)
y <- matrix(rnorm(24, mean = 0.2), nrow = 8, ncol = 3)
chen_qin_two_sample_test(x, y)

Chen–Song–Feng rank-based max white-noise test

Description

Implements the two rank procedures for which the primary paper supplies direct, scalable statistics and complete extreme-value calibrations: Spearman's rho and Kendall's tau. The paper does not establish rank-based sum or adaptive Fisher tests; those constructions in the book draft are not used here.

Usage

chen_song_feng_rank_white_noise_test(
  x,
  lag = 1L,
  measure = c("spearman", "kendall"),
  keep_lag = FALSE
)

Arguments

Value

An htest object. Exact ties are rejected because the published null variance and distribution-free calibration assume continuous margins.

References

Chen, D., Song, F. and Feng, L. Rank Based Tests for High Dimensional White Noise. Statistica Sinica 35, 1323–1347. doi:10.5705/ss.202022.0382.

Examples

t <- seq_len(18)
x <- cbind(sin(t), cos(t / 2), sin(t / 3 + 0.2))
chen_song_feng_rank_white_noise_test(x, lag = 2, measure = "spearman")

Chen–Zhang–Zhong high-dimensional covariance test

Description

Constructs the location-invariant unbiased estimators

T_{1n}=Y_{1n}-Y_{3n},\qquad T_{2n}=Y_{2n}-2Y_{4n}+Y_{5n}

of \operatorname{tr}(\Sigma) and \operatorname{tr}(\Sigma^2). For sphericity it uses U_n=pT_{2n}/T_{1n}^2-1; for identity it uses V_n=T_{2n}/p-2T_{1n}/p+1. In either case the primary null calibration is nU_n/2 or nV_n/2, respectively, against an upper standard-normal tail.

Usage

chen_zhang_zhong_covariance_test(x, null = c("sphericity", "identity"))

Arguments

Value

An hd_covariance_test object retaining all five U-statistic building blocks.

References

Chen, S. X., Zhang, L.-X. and Zhong, P.-S. (2010). Journal of the American Statistical Association, 105, 810–819. doi:10.1198/jasa.2010.tm09560.

Examples

set.seed(42)
x <- matrix(rnorm(48), nrow = 12, ncol = 4)
chen_zhang_zhong_covariance_test(x, null = "sphericity")

Cheng–Liu–Peng–Zhang–Zheng two-sample SSCM test

Description

Tests equality of two population spatial-sign covariance matrices (SSCMs), equivalently proportionality of the two scatter matrices under elliptical symmetry. With separately estimated spatial medians and centered signs, the feasible statistic is

T=p(A+B-2C),

where A and B average squared inner products over ordered within-group pairs and C averages them over all cross pairs.

Usage

cheng_sscm_equality_test(
  x,
  y,
  alpha = 0.05,
  tol = 1e-08,
  max_iter = 500L,
  zero_tol = 0,
  strict = TRUE
)

Arguments

Details

The implementation uses the complete finite-sample bias expression in Cheng et al. (2019), Equation (3.3). For each group it estimates the inverse radial moments through order three, evaluates the n^{-2} and n^{-3} terms separately, and subtracts p\widehat\delta. Under the null, Remark 3.1 estimates

q=p^{-1}\mathrm{tr}(\Lambda^2) =\frac{p}{n_1+n_2}\{n_1(A-\delta_1)+n_2(B-\delta_2)\}.

Consequently the standard-error denominator and its square are

\widehat v_0= \frac{2(n_1^{-1}+n_2^{-1})}{p+2}(p\widehat q),\qquad \widehat v_0^2=\{\widehat v_0\}^2.

The reported statistic is Z=(T-p\widehat\delta)/\widehat v_0 and the paper's rejection rule is the upper standard-normal tail.

A non-positive trace estimate, a zero fitted radius, or an unconverged spatial median invalidates calibration. No absolute value, variance floor, radius perturbation, or replacement center is used. With strict = FALSE the raw components are returned, but statistic, p.value, and the rejection decision are NA and diagnostics$calibrated is FALSE. The plug-in bias calibration has the additional regime stated in Remark 3.2 of the primary paper: under bounded eigenvalues it requires p=o(n_l^{3/2}) for both groups. This is stronger than merely using the theorem's p=O(n_l^2) expansion condition.

Value

An object of class c("hd_proportionality_test", "htest"). components contains all ordered-pair components, bias terms, scaled inverse radial moments, the trace estimate, the standard-error denominator and its squared variance, fitted centers, signs, and SSCMs.

References

Cheng, G., Liu, B., Peng, L., Zhang, B., and Zheng, S. (2019). Testing the equality of two high-dimensional spatial sign covariance matrices. Scandinavian Journal of Statistics, 46, 257–271. doi:10.1111/sjos.12350.

Examples

x <- matrix(c(-2, 0, 1, 2, 1, -1, 0, 2, -1, -2, 2, 0), 6, 2)
y <- matrix(c(-1, 1, 2, -2, 0, 3, 1, -1, 2, 0, -2, 1), 6, 2)
cheng_sscm_equality_test(x, y)

CHIME clustering for a two-component Gaussian mixture

Description

Fits the CHIME iteration of Cai, Ma, and Zhang for a two-component Gaussian mixture with a common covariance. At stage t, the sparse discriminant vector is the certified solution of

\min_\beta\;\tfrac12\beta^T\widehat\Sigma_t\beta- \beta^T(\widehat\mu_{1t}-\widehat\mu_{2t})+ \lambda_t\|\beta\|_1.

The covariance is never inverted. Coordinate descent uses the actual diagonal of the quadratic metric and must satisfy the exact lasso KKT equations before an iterate is accepted.

Usage

chime_clustering(
  x,
  lambda,
  initial,
  K = 2L,
  max_iter = length(lambda) - 1L,
  lasso_tol = 1e-08,
  lasso_max_iter = 10000L,
  support_tol = 0,
  psd_tol = sqrt(.Machine$double.eps)
)

Arguments

Details

This function deliberately differs from the authors' numerical script in three ways required for a reusable unsupervised method: it applies no hidden covariance ridge, it never uses true labels to choose a penalty, and it requires the complete penalty path (including the initial beta fit) from the caller. The method is only defined here for K = 2, matching the stated model and the available primary implementation.

Value

A chime_fit object containing labels, responsibilities, all parameter stages, the explicit lambda path, covariance matrices, and per-stage objective, convergence, and KKT certificates. Score ties use the printed weak inequality and are assigned to component 1.

References

Cai, T. T., Ma, J., and Zhang, L. (2019). CHIME: Clustering of high-dimensional Gaussian mixtures with EM algorithm and its optimality. Annals of Statistics, 47, 1234–1267. doi:10.1214/18-AOS1711.

Examples

x <- rbind(
  c(-2.2, -1.0), c(-1.8, -1.2), c(-2.0, -0.8), c(-1.7, -1.1),
  c( 1.8,  1.0), c( 2.2,  1.1), c( 2.0,  0.9), c( 1.7,  1.2)
)
initial <- list(
  omega = 0.5, mu1 = c(-1.5, -0.8), mu2 = c(1.5, 0.8),
  covariance = diag(2)
)
chime_clustering(x, lambda = c(0.2, 0.15, 0.1), initial = initial,
                 max_iter = 2)

Classical canonical correlation analysis

Description

Forms all three covariance blocks with one explicit divisor, constructs the symmetric inverse square roots of the two within-block covariance matrices, and applies an SVD to the whitened cross-covariance matrix. Both within-block matrices must be strictly positive definite and full rank at rank_tol. No ridge, pseudoinverse, or eigenvalue floor is used.

Usage

classical_cca(
  x,
  y,
  components = NULL,
  center = c("mean", "none"),
  covariance_divisor = c("n-1", "n"),
  rank_tol = sqrt(.Machine$double.eps)
)

Arguments

Value

A classical_cca_fit object with canonical coefficients, scores, structure loadings, correlations, covariance blocks, and whitened singular-vector projectors.

References

Hotelling, H. (1936). Relations between two sets of variates. Biometrika, 28, 321–377.

Examples

basis <- stats::poly(seq_len(8), degree = 4)
x <- basis[, 1:2, drop = FALSE]
y <- cbind(0.8 * basis[, 1] + 0.6 * basis[, 3],
           0.3 * basis[, 2] + sqrt(0.91) * basis[, 4])
colnames(x) <- c("x1", "x2")
colnames(y) <- c("y1", "y2")
fit <- classical_cca(x, y, components = 1)
fit$canonical.correlations

Classical CUSUM change-point test

Description

Computes the univariate CUSUM

C_n(k)=n^{-1/2}\{S_k-kS_n/n\}

and rejects for a large trimmed maximum of \{(k/n)(1-k/n)\}^{-\gamma}|C_n(k)|/\widehat\sigma. For gamma = 0, calibration = "asymptotic" uses the classical two-sided Brownian-bridge (Kolmogorov) limit. calibration = "grid-gaussian" simulates the method's finite scan-grid Gaussian bridge; this is null calibration, not replication of a paper simulation study.

Usage

classical_cusum_test(
  x,
  gamma = 0,
  trim = 1L,
  sigma = NULL,
  variance = c("difference", "sample"),
  calibration = c("asymptotic", "grid-gaussian"),
  calibration_draws = 4999L,
  alpha = 0.05,
  seed = NULL,
  strict = TRUE
)

Arguments

Details

If sigma is not supplied, the scale is either the sample standard deviation or Rice's adjacent-difference estimate. This software choice is returned explicitly because the book's classical display assumes a population \sigma and does not prescribe its estimator.

Value

An htest-compatible object with the complete CUSUM path.

References

Page (1954); Csorgo and Horvath (1997).

Examples

x <- c(-1, 0, 1, -1, 0, 3, 4, 3)
classical_cusum_test(x)

Classical plug-in linear discriminant classifier

Description

Estimates the common covariance with the pooled unbiased estimator and inserts it into the Gaussian LDA log-likelihood ratio. The covariance must be strictly positive definite. No generalized inverse or ridge is used.

Usage

classical_lda_classifier(x, y, prior = "equal", strict = TRUE)

Arguments

Value

An hd_classifier_fit.

References

Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis, 3rd ed. Wiley.

Examples

x <- rbind(c(2, 0), c(1, 1), c(1, -1),
           c(-2, 0), c(-1, 1), c(-1, -1))
classical_lda_classifier(x, rep(c("A", "B"), each = 3))

Classical principal component analysis with explicit covariance convention

Description

Computes PCA from the centered second-moment matrix with a user-visible divisor. The book's displayed estimator uses covariance_divisor = "n"; "n-1" gives the usual unbiased sample covariance. No variance floor or rank repair is applied. components = NULL returns all numerically usable positive-variance components.

Usage

classical_pca(
  x,
  components = NULL,
  center = c("mean", "none"),
  covariance_divisor = c("n-1", "n"),
  eigen_tol = sqrt(.Machine$double.eps)
)

Arguments

Details

Eigenvector signs are anchored by making the first coordinate attaining the largest absolute loading positive. Individual eigenvectors inside a repeated-eigenvalue block are not identified; the returned projector is the invariant object for such a block.

Value

A classical_pca_fit and hd_pca_fit object containing loadings, scores, covariance operator, eigenvalues, projector, reconstruction, and complete centering/divisor diagnostics.

References

Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis, 3rd ed. Wiley.

Examples

x <- rbind(c(2, 0), c(1, 1), c(-1, -1), c(-2, 0))
fit <- classical_pca(x, components = 1, covariance_divisor = "n")
fit$loadings

Classical plug-in quadratic discriminant classifier

Description

Uses classwise unbiased covariance matrices with divisors n1 - 1 and n2 - 1, then forms the canonical Gaussian QDA log-likelihood ratio. Both covariance matrices must be strictly positive definite. Ordinary positive determinants are used; no pseudo-determinant, absolute value, ridge, or generalized inverse is substituted.

Usage

classical_qda_classifier(x, y, prior = "equal", strict = TRUE)

Arguments

Value

An hd_classifier_fit.

References

Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis, 3rd ed. Wiley.

Examples

x <- rbind(c(2, 0), c(1, 1), c(1, -1), c(2, 2),
           c(-2, 0), c(-1, 2), c(-1, -2), c(-2, -1))
classical_qda_classifier(x, rep(c("A", "B"), each = 4))

Classical spatial sign and rank location tests

Description

Implements the fixed-dimensional spatial-sign test, the one-sample spatial signed-rank test, and the pooled two-sample spatial-rank test. All three use the zero-direction convention U(0)=0 and a Cholesky solve; a singular studentizing matrix is reported rather than replaced by a generalized inverse.

Usage

spatial_sign_test(x, mu = NULL, tol = sqrt(.Machine$double.eps), zero_tol = 0)

spatial_signed_rank_test(
  x,
  mu = NULL,
  tol = sqrt(.Machine$double.eps),
  zero_tol = 0
)

spatial_rank_test(x, y, tol = sqrt(.Machine$double.eps), zero_tol = 0)

Arguments

Details

For centered observations Z_i=X_i-\mu_0, the sign statistic uses

\bar U=n^{-1}\sum_i U(Z_i),\quad B_U=n^{-1}\sum_i U(Z_i)U(Z_i)^T,\quad Q_{sign}=n\bar U^T B_U^{-1}\bar U.

The signed ranks are

R_i=n^{-1}\sum_j U(Z_i+Z_j).

Since the first-order projection of their average has covariance 4B_R, where B_R=n^{-1}\sum_i R_iR_i^T, its statistic is

Q_{SR}=\frac{n}{4}\bar R^T B_R^{-1}\bar R.

The factor 1/4 is essential.

For two samples, every observation is ranked against the pooled sample:

R(Y_i)=N^{-1}\sum_j U(Y_i-Y_j).

If C=(N-1)^{-1}\sum_i R(Y_i)R(Y_i)^T, the sample covariance of the pooled ranks (whose average is exactly zero), the statistic is

Q_{2SR}=\frac{n_1n_2}{N}(\bar R_1-\bar R_2)^T C^{-1}(\bar R_1-\bar R_2).

These calibrations are asymptotic chi-squared laws for fixed dimension, not high-dimensional approximations. Euclidean spatial signs make the tests invariant to translations, common positive rescaling, and orthogonal transformations, but not to arbitrary nonspherical affine transformations.

Value

An object of classes hd_location_test and htest. statistic is the chi-squared statistic, parameter is its asymptotic degrees of freedom, variance is the studentizing second-moment matrix (or ⁠4 B_R⁠ for the signed-rank test), and components contains the directional mean and unscaled moment matrices. Numerical rank, reciprocal condition number, solver, zero counts, and calibration type are in diagnostics.

References

Mottonen, J. and Oja, H. (1995). Multivariate spatial sign and rank methods. Journal of Nonparametric Statistics, 5, 201–213.

Oja, H. (2010). Multivariate Nonparametric Methods with R. Springer.

Examples

set.seed(12)
x <- matrix(rnorm(80), 20, 4)
spatial_sign_test(x)
spatial_signed_rank_test(x)

y <- matrix(rnorm(96, 0.25), 24, 4)
spatial_rank_test(x, y)

Gaussian CLIME sparse precision estimator

Description

Solves the Cai–Liu–Luo CLIME program column by column,

\min_b\|b\|_1\quad\text{subject to}\quad \|S_n b-e_j\|_\infty\leq\lambda,

using a primal–dual algorithm. Every raw column must pass primal feasibility, dual feasibility, lasso stationarity, and relative duality-gap certificates. The final symmetric estimator retains the entry of smaller absolute value from each transposed pair; exact ties retain the row-column entry before mirroring, making the convention deterministic.

Usage

clime_precision(
  x,
  lambda,
  center = TRUE,
  divisor = c("n", "n-1"),
  solver_tol = 1e-07,
  solver_max_iter = 100000L,
  strict = TRUE
)

Arguments

Details

Primary CLIME does not guarantee that this post-LP symmetrization remains feasible, nor that the final matrix is positive definite. Both facts are reported and neither is repaired. Consequently valid certifies the raw column programs, not an invented SPD condition. SCIO and scaled-lasso precision estimation remain review-only because their additional programs are not specified by the short book display.

Value

A gaussian_precision_fit list. estimate is the primary smaller-absolute-value symmetrization of the certified raw columns.

References

Cai, T. T., Liu, W. and Luo, X. (2011). A constrained \ell_1 minimization approach to sparse precision matrix estimation. Journal of the American Statistical Association, 106, 594–607. doi:10.1198/jasa.2011.tm10155.

Examples

x <- rbind(c(-2, 0), c(-1, -1), c(1, 1), c(2, 0))
clime_precision(x, lambda = 0)

Feng–Zou–Wang–Zhu Composite T-squared two-sample test

Description

Implements the two-sample Composite T^2 test of Feng, Zou, Wang, and Zhu (2017). Observations are rows and variables are columns. Despite the wording at ch2_location.tex:778--784 in the current book draft, this is not a Behrens–Fisher test: the paper assumes that the two groups have the same unknown covariance matrix. Nor does the method combine existing test p-values or estimate correlations among component tests. Instead, it builds a block-diagonal approximation to the common pooled covariance and combines the resulting block Hotelling quadratic forms.

Usage

composite_t2_two_sample_test(x, y, block_size = 2L, selection = "paper_greedy")

Arguments

Details

For a fixed leave-out covariance, the paper's practical block rule starts with the pair having the largest absolute sample correlation and repeatedly adds the variable with the largest sum of absolute correlations to the current block. Variables already assigned to a block are removed and the procedure is repeated. selection = "paper_greedy" implements that rule, with lexicographic tie-breaking. The final block contains all remaining variables when fewer than block_size remain. The paper's notation and simulations use equal-size blocks; this deterministic remainder convention extends the displayed rule to dimensions not divisible by the block size.

Let \widehat{\Sigma}_{O_K,i_1,i_2,j_1,j_2}^{-1} be the block-diagonal inverse obtained after deleting observations i_1,i_2 from group 1 and j_1,j_2 from group 2, recomputing the unbiased pooled covariance and its correlation matrix, and rebuilding the blocks. The published statistic is

Q_n={1\over n_1n_2(n_1-1)(n_2-1)} \sum_{i_1\ne i_2}\sum_{j_1\ne j_2} (X_{1i_1}-X_{2j_1})^T \widehat{\Sigma}_{O_K,i_1,i_2,j_1,j_2}^{-1} (X_{1i_2}-X_{2j_2}).

The upper-tail normal calibration is

Z={Q_n\over {2(n_1^{-1}+n_2^{-1})^2 \widehat{\mathrm{tr}(\Lambda_K^2)}}^{1/2}}.

Following Theorem 2 and the paragraph immediately after it, the feasible trace estimate is the paper's one-sample leave-four-out estimator computed from group 1 only:

\widehat{\mathrm{tr}(\Lambda_K^2)}= {1\over 2P_{n_1}^4}\sum^* (X_{i_1}-X_{i_2})^T\widehat{\Sigma}_{O_K,\setminus4}^{-1} (X_{i_3}-X_{i_4}) (X_{i_1}-X_{i_4})^T\widehat{\Sigma}_{O_K,\setminus4}^{-1} (X_{i_3}-X_{i_2}).

Thus Q_n and its population scaling are symmetric in the two samples, but the paper's feasible finite-sample standardisation is label-asymmetric. Exchanging the samples can change the reported Z and p-value.

The first group must satisfy nrow(x) >= block_size + 5: after leaving out four rows, a block of size block_size then has enough residual degrees of freedom to be nonsingular. The second group must contain at least three rows for the pooled 2+2 leave-out covariance. These count conditions are necessary, not sufficient. Every marginal variance used for selection and every selected full or leave-out covariance block must be finite and strictly positive definite. Failure produces an error. No ridge, generalized inverse, absolute-value repair, or variance floor is applied.

The C++ kernel uses sufficient statistics and sums unordered pairs and quadruples as exact symmetry reductions of the displayed ordered formulas. Both samples are internally translated by a common anchor and divided by a common positive scale in each column. This is algebraically neutral and protects the scale-invariant method under very large or small units.

Value

An object of class c("hd_location_test", "htest"). The statistic is the standardised Z and the p-value is its upper standard-normal tail probability. components retains Q_n, the group-1 leave-four-out trace estimate, every calibration coefficient and ordered denominator, combination counts, and a full-sample diagnostic partition. diagnostics records the common-covariance assumption, dynamic block selection, first-sample calibration, conditioning summaries, and the absence of numerical repair.

References

Feng, L., Zou, C., Wang, Z., and Zhu, L. (2017). Composite T-squared test for high-dimensional data. Statistica Sinica, 27, 1419–1436. doi:10.5705/ss.202015.0199.

Examples

set.seed(2717)
x <- matrix(rnorm(21), 7, 3)
y <- matrix(rnorm(15, 0.2), 5, 3)
composite_t2_two_sample_test(x, y)

Construct a conditional-alpha sieve design from a supplied basis

Description

This helper performs only the algebra specified by the conditional-alpha papers. Rows are time observations. If center_alpha = TRUE, its first block is the column-centered supplied basis; otherwise it is the original basis. The remaining blocks are, in factor order, the rowwise products f_{jt}B(t/T). It does not choose spline order, knots, basis dimension, or a BIC rule.

Usage

conditional_alpha_sieve_design(
  basis,
  factors = NULL,
  center_alpha = TRUE,
  alpha_contrast = NULL
)

Arguments

Details

A normalized B-spline basis usually contains the constant function, so all of its centered columns are linearly dependent. The papers write an ordinary inverse despite this identity. This function never silently drops a column: pass an explicit full-column-rank alpha_contrast, or supply a reduced basis. The contrast changes only the coordinates, not the spanned centered-alpha space, when it has the intended range.

Value

A full-column-rank design matrix carrying its basis/factor metadata.

References

Ma, S., Lan, W., Su, L. and Tsai, C.-L. (2020). Testing alphas in conditional time-varying factor models with high-dimensional assets. Journal of Business & Economic Statistics, 38, 214–227. doi:10.1080/07350015.2018.1482758.

Examples

tt <- seq(0, 1, length.out = 12)
basis <- cbind(tt, tt^2)
factors <- cbind(market = sin(seq_len(12)))
conditional_alpha_sieve_design(basis, factors)

Fit the restricted conditional-alpha sieve regression

Description

Fits, independently for every asset, the no-intercept null regression on a supplied full-rank sieve design. The residualized intercept h=M_Z1_T is retained because it enters every feasible calibration. The function uses one common response scale and independent design-column scales for numerical stability; both leave the fitted residual subspace and every reported test unchanged. It never uses a generalized inverse or ridge.

Usage

conditional_alpha_sieve_fit(returns, design)

Arguments

Value

A reusable conditional_alpha_sieve_fit object containing the restricted residuals, h, design metadata, and projection diagnostics.

References

Ma, S., Lan, W., Su, L. and Tsai, C.-L. (2020). Testing alphas in conditional time-varying factor models with high-dimensional assets. Journal of Business & Economic Statistics, 38, 214–227. doi:10.1080/07350015.2018.1482758.

Examples

tt <- seq(0, 1, length.out = 14)
z <- conditional_alpha_sieve_design(cbind(tt, tt^2))
y <- cbind(sin(1:14), cos(1:14), sin(1:14 / 2))
conditional_alpha_sieve_fit(y, z)

Fixed-dimensional conditional-factor Wald benchmark

Description

Evaluates the Chapter 4 supplied-estimate benchmark

W=T\hat\delta^T\hat\Omega_\delta^{-1}\hat\delta,

with a chi-squared reference distribution having length(delta) degrees of freedom. This function does not estimate spline or kernel nuisances.

Usage

conditional_factor_wald_test(delta, covariance, sample_size)

Arguments

Value

A strict Cholesky-based upper-tail Wald htest object.

References

Li, D. and Yang, L. (2011). Nonparametric tests of conditional factor models. Ang, A. and Kristensen, D. (2012). Testing conditional factor models.

Examples

conditional_factor_wald_test(c(0.1, -0.2), diag(c(1, 2)), 40)

Gaussian finite-power moments for the analytical aSPU test

Description

Gaussian finite-power moments for the analytical aSPU test

Usage

cpp_aspu_power_moments(scores, covariance, powers)

Arguments

Value

Internal list of observed sums and Gaussian moment matrices.

Construct standardised contrasts for the analytical aSPU test

Description

Construct standardised contrasts for the analytical aSPU test

Usage

cpp_aspu_standardize(
  x,
  y,
  correlation_source,
  bandwidth1,
  bandwidth2,
  supplied_correlation,
  supplied_standard_errors
)

Arguments

Value

Internal list of contrasts and covariance diagnostics.

Original Bai–Saranadasa two-sample statistic kernel

Description

Original Bai–Saranadasa two-sample statistic kernel

Usage

cpp_bai_saranadasa_two_sample(x, y)

Arguments

Value

Internal list of statistic components.

Spatial geometry and unweighted diagonal HR update

Description

Spatial geometry and unweighted diagonal HR update

Usage

cpp_ch2_generic_weighted_geometry(x, location, diagonal, zero_tol)

Arguments

Value

Directions, radii, and the literal unweighted HR scale update.

Strict initial values for the generic weighted HR recursion

Description

Strict initial values for the generic weighted HR recursion

Usage

cpp_ch2_generic_weighted_initial(x)

Arguments

Value

Column means and unbiased marginal variances.

Generic oracle weighted-sign quadratic U-statistic

Description

Generic oracle weighted-sign quadratic U-statistic

Usage

cpp_ch2_generic_weighted_quadratic(directions, weights)

Arguments

Value

The literal pairwise score and empirical weight moment.

Literal generic weighted HR location update

Description

Literal generic weighted HR location update

Usage

cpp_ch2_generic_weighted_step(location, diagonal, directions, radii, weights)

Arguments

Value

The weighted numerator, denominator, and next location.

Robust conditional spatial-sign sum ingredients

Description

Robust conditional spatial-sign sum ingredients

Usage

cpp_ch4_afc_css_components(residuals, trace_residuals, h)

Arguments

Value

CSS numerator, trace estimator, signs, and diagnostics.

Light-tail conditional-alpha sum and maximum ingredients

Description

Light-tail conditional-alpha sum and maximum ingredients

Usage

cpp_ch4_afc_light_components(
  residuals,
  h,
  factor_count,
  design_columns,
  compute_max
)

Arguments

Value

Exact feasible sum, trace, variance, and maximum ingredients.

Strict least-squares projection for conditional-alpha procedures

Description

Strict least-squares projection for conditional-alpha procedures

Usage

cpp_ch4_afc_project(y, design)

Arguments

Value

Restricted residuals, residualized intercept, coefficients, and exact numerical diagnostics.

Spatial Kendall matrix for robust latent-factor extraction

Description

Spatial Kendall matrix for robust latent-factor extraction

Usage

cpp_ch4_afc_spatial_kendall(x)

Arguments

Value

Spatial Kendall matrix and the number of zero pair differences.

Core OLS quantities for unconditional factor-pricing alpha tests

Description

Core OLS quantities for unconditional factor-pricing alpha tests

Usage

cpp_ch4_alpha_ols(y, factors)

Arguments

Value

A list of OLS quantities in the scaled coordinates.

Standardized residual radii for the Chapter 4 completion methods

Description

Standardized residual radii for the Chapter 4 completion methods

Usage

cpp_ch4_completion_standardized_radii(residuals, diagonal)

Arguments

Value

A vector of stable Euclidean radii.

Exact degenerate rank-U core for vector independence

Description

Exact degenerate rank-U core for vector independence

Usage

cpp_ch4_completion_vector_u_core(
  x,
  y,
  measure,
  permutations,
  max_kernel_evaluations,
  keep_estimates,
  keep_permutation
)

Arguments

Value

Exact observed and intrinsic-permutation components.

Weighted spatial-sign alpha quadratic form

Description

Weighted spatial-sign alpha quadratic form

Usage

cpp_ch4_completion_weighted_alpha_q(directions, h, weights)

Arguments

Value

The quadratic form and empirical second weight moment.

Liu–Feng–Ma spatial-sign alpha core

Description

Liu–Feng–Ma spatial-sign alpha core

Usage

cpp_ch4_lfm_spatial_sign_core(
  y,
  factors,
  tol,
  max_iter,
  zero_tol,
  compute_trace
)

Arguments

Value

Components of the primary statistic.

Euclidean assignment kernel for Chapter 7 spatial clustering

Description

Euclidean assignment kernel for Chapter 7 spatial clustering

Usage

cpp_ch7sc_assign_euclidean(x, centers, active, keep_matrix)

Arguments

Value

Internal deterministic assignment and distance diagnostics.

SSCM metric assignment kernel for Chapter 7

Description

SSCM metric assignment kernel for Chapter 7

Usage

cpp_ch7sc_assign_metric(x, centers, inverse_metric, keep_matrix)

Arguments

Value

Internal deterministic assignment and distance diagnostics.

Across-center feature scores for Sparse–SM

Description

Across-center feature scores for Sparse–SM

Usage

cpp_ch7sc_feature_scores(centers)

Arguments

Value

Internal center average and hard-screening scores.

Active-subspace geometry for Sparse–SM selectors

Description

Active-subspace geometry for Sparse–SM selectors

Usage

cpp_ch7sc_geometry(x, centers, labels, overall, active)

Arguments

Value

Internal between- and within-spatial-median geometry.

SSCM and exact inverse for Chapter 7 spatial clustering

Description

SSCM and exact inverse for Chapter 7 spatial clustering

Usage

cpp_ch7sc_sscm_metric(x, centers, labels, lambda, zero_tol, keep_signs)

Arguments

Value

Internal SSCM, inverse, and finite-sample certificates.

Chen–Qin two-sample statistic kernel

Description

Chen–Qin two-sample statistic kernel

Usage

cpp_chen_qin_two_sample(x, y)

Arguments

Value

Internal list of statistic components.

Adaptive-threshold precision estimator for the feasible CLX test

Description

The pooled covariance and the variance estimates used in its entrywise thresholds both have denominator N = n1 + n2. The threshold is delta * sqrt(theta_ij * log(p) / N). A relative eigenvalue floor is applied only when needed to invert the thresholded covariance matrix.

Usage

cpp_clx_adaptive_precision(x, y, delta, relative_eigen_floor)

Arguments

Value

Internal list containing the precision estimate and diagnostics.

Oracle or feasible Cai–Liu–Xia statistic kernel

Description

Oracle or feasible Cai–Liu–Xia statistic kernel

Usage

cpp_clx_two_sample(x, y, precision, oracle)

Arguments

Value

Internal list of statistic components.

Feng–Zou–Wang–Zhu Composite T-squared two-sample kernel

Description

Feng–Zou–Wang–Zhu Composite T-squared two-sample kernel

Usage

cpp_composite_t2_two_sample(x, y, block_size = 2L, selection = "paper_greedy")

Arguments

Value

Internal list of statistic, trace, and partition diagnostics.

Feng–Sun one-sample spatial-sign statistic kernel

Description

Feng–Sun one-sample spatial-sign statistic kernel

Usage

cpp_feng_sun_one_sample(x, mu, tolerance, max_iterations)

Arguments

Value

Internal list of feasible statistic components and leave-out fits.

Feng–Wang PDQ spatial-sign kernel

Description

Feng–Wang PDQ spatial-sign kernel

Usage

cpp_feng_wang_pdq_two_sample(
  x,
  y,
  quantile_prob,
  level,
  B,
  seed,
  keep_bootstrap,
  tolerance,
  max_iterations
)

Arguments

Value

A list of the statistic, nuisance fits, bootstrap calibration, and numerical diagnostics.

Feng–Zou–Wang two-sample multivariate-sign kernel

Description

Feng–Zou–Wang two-sample multivariate-sign kernel

Usage

cpp_feng_zou_wang_two_sample_sign(x, y, tolerance, max_iterations)

Arguments

Value

Internal list containing the feasible statistic and diagnostics.

Feng–Zou–Wang–Zhu Behrens–Fisher statistic kernel

Description

Feng–Zou–Wang–Zhu Behrens–Fisher statistic kernel

Usage

cpp_fzwz_bf_two_sample(x, y)

Arguments

Value

Internal list of statistic components.

Feng–Liu–Ma one-sample inverse norm sign statistic kernel

Description

Feng–Liu–Ma one-sample inverse norm sign statistic kernel

Usage

cpp_inst_one_sample(x, mu, tol, max_iter, zero_tol)

Arguments

Value

Internal list of primary statistic, direct feasible variance, cross-fitted nuisance diagnostics, and iteration diagnostics.

Li–Wang–Zou simpler spatial-sign two-sample kernel

Description

Li–Wang–Zou simpler spatial-sign two-sample kernel

Usage

cpp_li_wang_zou_two_sample_sign(x, y, tol, max_iter, zero_tol)

Arguments

Value

Internal list containing the feasible SST and diagnostics.

Park–Ayyala one-sample statistic kernel

Description

Park–Ayyala one-sample statistic kernel

Usage

cpp_park_ayyala_one_sample(x, mu)

Arguments

Value

Internal list of statistic components.

Full-sample scaled spatial median kernel

Description

Full-sample scaled spatial median kernel

Usage

cpp_scaled_spatial_median(x, target, center_on_target, tol, max_iter, zero_tol)

Arguments

Value

Internal fit, max statistic ingredients, and diagnostics.

Srivastava–Katayama–Kano statistic kernel

Description

Srivastava–Katayama–Kano statistic kernel

Usage

cpp_skk_two_sample(x, y)

Arguments

Value

Internal list of statistic components.

Srivastava–Du one-sample statistic kernel

Description

Srivastava–Du one-sample statistic kernel

Usage

cpp_srivastava_du_one_sample(x, mu)

Arguments

Value

Internal list of statistic components.

Huang–Liu–Zhou–Feng two-sample inverse norm sign kernel

Description

Huang–Liu–Zhou–Feng two-sample inverse norm sign kernel

Usage

cpp_tinst_two_sample(x, y, tol, max_iter, zero_tol)

Arguments

Value

Internal list of tINST statistic and convergence components.

Wang–Peng–Li one-sample spatial-sign statistic kernel

Description

Wang–Peng–Li one-sample spatial-sign statistic kernel

Usage

cpp_wang_peng_li_one_sample(x, mu)

Arguments

Value

Internal list of statistic and cross-validation components.

Wang–Xu approximate-randomization kernel

Description

Wang–Xu approximate-randomization kernel

Usage

cpp_wang_xu_approx_randomization(
  x,
  y,
  calibration,
  B,
  seed,
  max_exact,
  keep_randomized
)

Arguments

Value

Internal observed, pseudo-sample, reference, and diagnostic fields.

Weighted scaled spatial median and diagonal scale kernel

Description

Weighted scaled spatial median and diagonal scale kernel

Usage

cpp_weighted_scaled_spatial_median(x, m, tol, max_iter, zero_tol)

Arguments

Value

Internal list containing the fitted location, diagonal scale, radial quantities, and iteration diagnostics.

Yan–Zhao–Feng weighted max statistic kernel

Description

Yan–Zhao–Feng weighted max statistic kernel

Usage

cpp_yzf_weighted_max(x, mu, m, tol, max_iter, zero_tol)

Arguments

Value

Internal list containing the weighted max statistic and fit.

Yan–Zhao–Feng weighted max-sum statistic kernel

Description

Yan–Zhao–Feng weighted max-sum statistic kernel

Usage

cpp_yzf_weighted_maxsum(x, mu, m, tol, max_iter, zero_tol)

Arguments

Value

Internal list containing weighted max and feasible weighted sum components and all fit diagnostics.

One-sample Zhang–Feng marginal signed-rank scores

Description

One-sample Zhang–Feng marginal signed-rank scores

Usage

cpp_zhang_feng_one_sample_scores(x, mu)

Arguments

Value

Internal list of marginal scores and exact untied null moments.

Ouyang Parzen long-run variance estimate for rank-square scores

Description

Ouyang Parzen long-run variance estimate for rank-square scores

Usage

cpp_zhang_feng_parzen_tau(score, lag)

Arguments

Value

Internal list of sample autocovariances, weights, and tau squared.

Two-sample Zhang–Feng marginal WMW scores

Description

Two-sample Zhang–Feng marginal WMW scores

Usage

cpp_zhang_feng_two_sample_scores(x, y)

Arguments

Value

Internal list of marginal scores and exact untied null moments.

Zhang–Zhou–Guo one-sample normal-reference kernel

Description

Zhang–Zhou–Guo one-sample normal-reference kernel

Usage

cpp_zhang_zhou_guo_one_sample(residual)

Arguments

Value

Internal list of the centred U-statistic, unbiased trace estimates, cumulants, and shifted-scaled chi-square parameters.

Direct sparse discriminant analysis

Description

Fits the Mai–Zou–Yuan DSDA lasso with the exact response coding -n/n_1 for class 1 and +n/n_2 for class 2 and objective

n^{-1}\sum_i(y_i-a-x_i^\top\beta)^2+ \lambda\|\beta\|_1.

The regression coefficient points toward class 2. The package reverses it only when constructing its documented score, so positive always means class 1. Following the primary equal-prior rule, classification uses the midpoint intercept rather than the penalized regression intercept. The coordinate-descent solution must pass both KKT and Fenchel dual-gap checks.

Usage

dsda_classifier(
  x,
  y,
  lambda,
  prior = c(0.5, 0.5),
  solver_tol = 1e-07,
  solver_max_iter = 100000L,
  strict = TRUE,
  tie = c("class1", "class2")
)

Arguments

Value

An object of class hd_classifier_fit.

References

Mai, Q., Zou, H. and Yuan, M. (2012). A direct approach to sparse discriminant analysis in ultra-high dimensions. Biometrika, 99, 29–42. doi:10.1093/biomet/asr066.

Examples

x <- rbind(
  c(2, 1), c(1, 2), c(2, 2), c(3, 1),
  c(-2, -1), c(-1, -2), c(-2, -2), c(-3, -1)
)
y <- factor(rep(c("first", "second"), each = 4))
dsda_classifier(x, y, lambda = 0.2)

EC2 sparse covariance estimation with an eigenvalue constraint

Description

Implements the convex off-diagonal-lasso branch of Liu, Wang and Zhao's EC2 estimator. The primary method first forms the empirical correlation matrix R, solves

\min_{C:\,\operatorname{diag}(C)=1,\,\lambda_{\min}(C)\geq\tau} \frac12\|R-C\|_F^2+\lambda\sum_{i\ne j}|c_{ij}|,

and restores marginal empirical standard deviations. This differs from the book's shortened direct-covariance display. The default covariance divisor is n; selecting "n-1" is explicit and changes only the marginal covariance scale, not the sample correlation.

Usage

ec2_covariance(
  x,
  lambda,
  tau,
  penalty = "l1",
  center = TRUE,
  divisor = c("n", "n-1"),
  rho = 1,
  solver_tol = 1e-07,
  solver_max_iter = 50000L,
  strict = TRUE
)

Arguments

Details

The ISP/ADMM solver uses the paper's soft-threshold and spectral-projection updates. A fit is valid only when the equality, fixed-point, off-diagonal subgradient, spectral-dual, complementarity, diagonal, and eigenvalue certificates all pass. The method-defining tau projection is not a numerical repair. Adaptive EC2 and MC+ EC2 require additional weight or shape choices and remain review-only; this function never guesses them.

Value

An ec2_covariance_fit list. estimate is the covariance matrix only for a certified fit; correlation is its EC2 correlation estimate.

References

Liu, H., Wang, L. and Zhao, T. (2014). Sparse covariance matrix estimation with eigenvalue constraints. Journal of Computational and Graphical Statistics, 23, 439–459. doi:10.1080/10618600.2013.782818.

Examples

x <- rbind(c(-2, -1), c(-1, 0), c(1, 0), c(2, 1))
ec2_covariance(x, lambda = 0.2, tau = 0.1)

Select the number of spikes in an elliptical factor model

Description

Computes the trace-p spatial-sign pilot centered at the sample spatial median and applies the eigenvalue-ratio (ER) or growth-ratio (GR) selector of Xu, Ma, Wang and Feng. For ordered eigenvalues \lambda_j,

\widehat m_{ER}=\arg\max_{1\le j\le M} \lambda_j/\lambda_{j+1}.

The GR criterion uses V_j=\sum_{l=j+1}^{\min(n,p)-1}\lambda_l exactly as in the primary paper. max_factors is deliberately supplied: the paper treats M as a predetermined upper bound and its numerical choice is not a universal tuning rule.

Usage

elliptical_factor_number(
  x,
  max_factors,
  method = c("er", "gr"),
  tol = 1e-08,
  max_iter = 500L,
  zero_tol = 0,
  strict = TRUE
)

Arguments

Value

An elliptical_factor_number object containing the selected count, all pilot eigenvalues, criterion values, and spatial-median diagnostics.

References

Xu, X., Ma, H., Wang, H. and Feng, L. (2025). High dimensional matrix estimation through elliptical factor models. arXiv:2512.19325. https://arxiv.org/abs/2512.19325.

Examples

x <- rbind(c(-3, -2, 0), c(-2, -1, 1), c(-1, 0, -1),
           c(1, 0, 1), c(2, 1, -1), c(3, 2, 0))
elliptical_factor_number(x, max_factors = 1)

Sparse precision estimation for an elliptical factor fit

Description

Applies the primary CLIME or GLASSO problem to the raw idiosyncratic complement of a spatial_sign_poet() or poet_tme() fit, and reconstructs the full inverse scatter through the Woodbury identity,

V=V_u-V_u\Gamma(\Lambda^{-1}+\Gamma^T V_u\Gamma)^{-1} \Gamma^T V_u.

CLIME must pass per-column primal, dual, stationarity, gap, and final symmetrized-feasibility certificates. GLASSO must remain positive definite and satisfy its full elementwise-\ell_1 KKT equation. The tuning parameter is explicit because the theoretical constant multiplying the rate is not an observable default.

Usage

elliptical_factor_precision(
  fit,
  lambda,
  method = c("sclime", "sglasso"),
  solver_tol = 1e-07,
  solver_max_iter = 100000L,
  initial_step = 1,
  max_backtracking = 100L,
  strict = TRUE
)

Arguments

Value

An elliptical_factor_precision_fit containing the idiosyncratic and full precision estimates and complete solver certificates.

References

Xu, X., Ma, H., Wang, H. and Feng, L. (2025). arXiv:2512.19325.

Examples

x <- rbind(c(-3, -2), c(-2, -1), c(-1, 1),
           c(1, -1), c(2, 1), c(3, 2))
f <- spatial_sign_poet(x, factors = 0, threshold = 0.2)
elliptical_factor_precision(f, lambda = 0.5, method = "sglasso")

Exact oracle classifier for two elliptical populations

Description

Computes the exact generator-aware score

\log\pi_1-\tfrac12\log|\Lambda_1|+\log g_1(d_1) -\log\pi_2+\tfrac12\log|\Lambda_2|-\log g_2(d_2),

where d_k=(z-\mu_k)'\Lambda_k^{-1}(z-\mu_k). The log-generator functions must be vectorized. They may return -Inf for a zero density, but not NA, NaN, or Inf.

Usage

elliptical_oracle_classifier(
  location1,
  location2,
  shape1,
  shape2,
  log_generator1,
  log_generator2 = log_generator1,
  prior = "equal",
  levels = c("class1", "class2"),
  feature_names = NULL
)

Arguments

Details

A common generator and common shape do not generally reduce this rule to a midpoint linear rule under unequal priors. That shortcut is exact for equal priors, and for the exponential radial generator underlying Gaussian LDA.

Value

A generator-aware hd_classifier_fit.

References

Fang, K.-T. and Anderson, T. W. (1990). Statistical Inference in Elliptically Contoured and Related Distributions. Allerton Press. Wakaki, H. (1994). Discriminant analysis under elliptical populations. Hiroshima Mathematical Journal, 24, 257–298.

Examples

normal_log_generator <- function(d) -d / 2
fit <- elliptical_oracle_classifier(
  c(1, 0), c(-1, 0), diag(2), diag(2), normal_log_generator
)
predict(fit, rbind(c(1, 0), c(-1, 0)))

Elliptical regularized Hotelling test with Cauchy aggregation

Description

Computes the feasible ERHT statistic over a fixed deterministic ridge grid and combines its upper-tail normal p-values with the analytic Cauchy rule of Feng, Zhou, and Wang (2026). Equal weights are used by default. The default grid \{0.1,0.2,\ldots,1\} is the grid used in the paper's numerical implementation; this function does not reproduce any simulation design from that paper.

Usage

elliptical_regularized_hotelling_cauchy_test(
  x,
  mu = NULL,
  rho = seq(0.1, 1, by = 0.1),
  weights = NULL,
  alpha = 0.05,
  tol = 1e-08,
  max_iter = 1000L,
  strict = TRUE,
  keep_companion = FALSE
)

Arguments

Details

For positive weights \varpi_k summing to one,

T_{CC}=\sum_k\varpi_k\tan[\pi\{1/2-p_k\}],\qquad p_{CC}=1/2-\pi^{-1}\arctan(T_{CC}).

The computation uses signed-log cotangents and an angle representation to avoid overflow and cancellation in extreme tails. The article proves dependence-robust small-tail validity for this analytic p-value; except in special dependence structures it does not claim exact fixed-level calibration at an ordinary level such as 0.05.

Value

An object of class c("hd_location_test", "htest"). Its p-value is the analytic Cauchy combination. The marginal table contains every raw ERHT statistic, feasible center and variance, standardized statistic, p-value, log-p-value, and fixed-ridge decision.

References

Feng, L., Zhou, L., and Wang, X. (2026). Elliptical regularized Hotelling testing for high dimensional data. arXiv:2606.25942. doi:10.48550/arXiv.2606.25942.

Examples

set.seed(2)
x <- matrix(rnorm(160), 40, 4)
elliptical_regularized_hotelling_cauchy_test(
  x, rho = c(0.2, 0.5, 1)
)

Elliptical regularized Hotelling fixed-ridge test

Description

Tests a high-dimensional location vector under elliptical symmetry and pervasive dependence using the fixed-ridge ERHT statistic of Feng, Zhou, and Wang (2026). Observations are rows and variables are columns.

Usage

elliptical_regularized_hotelling_test(
  x,
  mu = NULL,
  rho = 0.5,
  alpha = 0.05,
  tol = 1e-08,
  max_iter = 1000L,
  strict = TRUE,
  keep_companion = FALSE
)

Arguments

Details

Let \widehat\theta be the sample spatial median, \widehat Y_i=\sqrt p\,U(X_i-\widehat\theta), and

\widehat R_n=n^{-1}\sum_i\widehat Y_i\widehat Y_i^{\mathsf T}.

For \rho>0, the raw statistic is

T_n(\rho)=n(\widehat\theta-\theta_0)^{\mathsf T} (\widehat R_n+\rho I)^{-1}(\widehat\theta-\theta_0).

The function implements the paper's direct feasible companion-matrix centering n\widehat\mu_n(\rho) and variance n\widehat\sigma_{D,n}^2(\rho), returning the upper-tail normal calibration of

Z_n(\rho)=\{T_n(\rho)-n\widehat\mu_n(\rho)\}/ \{n\widehat\sigma_{D,n}^2(\rho)\}^{1/2}.

The C++ kernel uses one economy SVD of the sign matrix. It evaluates the ridge quadratic in row and null spaces and switches between the companion matrix A and its complement I-A for the feasible calibration. Thus it never forms a persistent p\times p inverse or subtracts nearly equal Woodbury, diagonal-weight, or companion terms. A zero fitted residual makes the required inverse distance undefined and is reported as an error; no observation is omitted or perturbed and no additional ridge, pseudoinverse, absolute-value repair, or variance floor is used.

The paper's numerical section mentions a separate Bartlett center correction with a=8n/p, but neither the article nor its public arXiv source specifies the formula that maps a into the feasible center. This function therefore implements the fully stated equations only and deliberately does not guess that simulation-program modification.

Value

An object of class c("hd_location_test", "htest"). The reported statistic is Z_n(\rho); raw.statistic is T_n(\rho). Components include all feasible functionals, raw and internally scaled values, inverse distances, spatial-median diagnostics, and companion eigenvalues. Raw dimensional fields can underflow to zero or overflow to infinity when their physical units exceed the double-precision range; the separately normalized fields in components$scaled remain the numerical-audit representation used to compute the reported statistic.

References

Feng, L., Zhou, L., and Wang, X. (2026). Elliptical regularized Hotelling testing for high dimensional data. arXiv:2606.25942. doi:10.48550/arXiv.2606.25942.

Examples

set.seed(1)
x <- matrix(rnorm(120), 30, 4)
elliptical_regularized_hotelling_test(x, rho = 0.5)

Elliptical regularized Hotelling change-point scan

Description

Implements the single-change path or the discretized adjacent-triple scan of Song, Wen, and Feng (2026). For adjacent segments with sizes n1, n2, the raw statistic is

V_\rho^{raw}=N\widehat\Delta^T (\widehat R+\rho I)^{-1}\widehat\Delta, \qquad N=n_1n_2/(n_1+n_2).

The common pool is the full sample. Its centered spatial signs are columns of Y; with A=Y^TQY/n, the primary studentization is

\widehat\kappa=\sum_i\beta_i^2A_{ii},\qquad \widehat\sigma^2=2n\sum_{i\ne j}\beta_i^2\beta_j^2A_{ij}^2,

Z_\rho=(V_\rho^{raw}-n\widehat\kappa)/ \{n\widehat\sigma^2\}^{1/2}.

The ordered off-diagonal variance is accumulated directly. It is never obtained by flooring a cancellation-prone difference.

Usage

erht_change_point_test(
  x,
  ridge,
  scan = c("single", "multiple"),
  epsilon = 0.1,
  calibration = c("gaussian-supremum", "time-permutation", "none"),
  calibration_draws = 4999L,
  alpha = 0.05,
  seed = NULL,
  tol = 1e-08,
  max_iter = 1000L,
  zero_tol = 0,
  cauchy_weights = NULL,
  keep_calibration = FALSE,
  strict = TRUE
)

Arguments

Details

ridge contains the actual positive \rho values in the primary formula. The public ERHTCP reproduction code instead accepts empirical ratios rho / (p/n); this function does not silently make that conversion. Users who want that convention can supply ridge_ratio * p / n explicitly, and all resolved ridge values are returned.

calibration = "gaussian-supremum" simulates the primary marginal Gaussian-process limit on the exact candidate grid. "time-permutation" is the paper's practical exchangeability calibration and preserves each multivariate row. It is invalid under unaddressed serial dependence. With multiple ridge values, the returned Cauchy transform combines the marginal p-values, but calibration.exact is FALSE: the paper explicitly distinguishes this analytic rule from exact calibration by a joint-limit quantile involving the unknown cross-ridge correlation r_E.

The primary inverse-distance average is exactly mean(sqrt(p) / distance), with a coincident observation assigned weight zero. The optional (p - 1) / p factor used by the public reproduction repository is not part of the displayed primary statistic and is not used.

Value

An htest object with every local raw/center/variance/Z component.

References

Song, Wen and Feng (2026), arXiv:2607.22162.

Examples

x <- rbind(
  c(-2, 0), c(2, 0), c(0, -2), c(0, 2),
  c(-1, -1), c(1, 1), c(-1, 1), c(1, -1),
  c(-2, 1), c(2, -1), c(-1, 2), c(1, -2)
)
erht_change_point_test(
  x, ridge = 0.5, epsilon = 0.25,
  calibration_draws = 99, seed = 1
)

Wild binary segmentation with the ERHT local score

Description

Implements the primary WBS–ERHT algorithm of Song, Wen, and Feng (2026). On an interval I, every (or, when triple_step > 1, an explicitly thinned) adjacent triple compares two segments of length at least ceiling(epsilon * length(I)); the common scatter pool is the whole interval. Scores are maximized over the supplied actual ridge grid.

Usage

erht_wbs(
  x,
  ridge,
  threshold,
  min_interval,
  refinement_radius,
  deletion_radius,
  intervals = NULL,
  M = NULL,
  epsilon = 0.1,
  triple_step = 1L,
  max_changes = Inf,
  seed = NULL,
  tol = 1e-08,
  max_iter = 1000L,
  zero_tol = 0,
  strict = TRUE
)

Arguments

Details

The recursion uses the paper's narrowest-over-threshold rule: among sampled intervals contained in the current segment, plus the current segment itself, retain those with score strictly greater than threshold; choose the narrowest, break a length tie by the larger score, refine within the requested radius, delete the requested neighborhood, and recurse.

The primary theory does not choose the ridge grid, WBS threshold, number of intervals, minimum length, refinement radius, or deletion radius. Hence this API requires those quantities (or explicit intervals) rather than inventing automatic tuning constants. Random WBS intervals are generated independently of the observations from uniformly sampled unordered endpoint pairs. Scanning/WBS is intrinsic to the estimator and is not replication of a paper simulation study.

triple_step = 1 is the literal exhaustive adjacent-triple collection. Larger values are an explicit computational approximation and are recorded. Any invalid spatial median or local variance invalidates the fit; failed triples are never silently removed.

Value

A change-point fit with sorted estimates and selection history.

References

Song, Wen and Feng (2026), arXiv:2607.22162, Algorithm 1.

Examples

x <- rbind(
  c(-2, 0), c(2, 0), c(0, -2), c(0, 2),
  c(-1, -1), c(1, 1), c(-1, 1), c(1, -1),
  c(-2, 1), c(2, -1), c(-1, 2), c(1, -2)
)
erht_wbs(
  x, ridge = 0.5, threshold = 10,
  min_interval = 8, refinement_radius = 3, deletion_radius = 1,
  intervals = matrix(c(1, 12), ncol = 2), epsilon = 0.25
)

Features annealed independence rule (FAIR)

Description

Implements Fan and Fan's FAIR classifier. Features are ranked by the absolute Welch two-sample statistic

T_j=(\bar X_j-\bar Y_j)/ \sqrt{S_{1j}^2/n_1+S_{2j}^2/n_2},

while classification uses the marginal variance (S_{1j}^2+S_{2j}^2)/2. Feature-count selection by "paper" maximizes equation (4.3), including the largest eigenvalue of every truncated pooled within-class correlation matrix. Ties in absolute t-statistics are broken by the original feature index, and a tie in the criterion selects the smallest feature count.

Usage

fair_classifier(
  x,
  y,
  selection = c("paper", "m", "threshold"),
  m = NULL,
  t_threshold = NULL,
  prior = "equal",
  zero_variance = c("error", "drop"),
  strict = TRUE
)

Arguments

Details

FAIR is implemented under its equal-prior common-covariance contract. An exactly zero marginal variance is never floored: it is either an error or is dropped explicitly.

Value

An hd_classifier_fit.

References

Fan, J. and Fan, Y. (2008). High-dimensional classification using features annealed independence rules. Annals of Statistics, 36, 2605–2637. doi:10.1214/07-AOS504.

Examples

x <- rbind(c(3, 1, 0), c(2, 2, 1), c(4, 0, -1),
           c(-3, -1, 0), c(-2, -2, -1), c(-4, 0, 1))
fair_classifier(x, rep(c("A", "B"), each = 3), selection = "m", m = 1)

Fantope projection and selection sparse PCA

Description

Solves

\max_H \langle S,H\rangle-\tau\lVert H\rVert_{1,1},\quad 0\preceq H\preceq I,\quad \mathrm{tr}(H)=r

by a certified two-block ADMM. The H update is the Euclidean Fantope projection and the split update is entrywise soft thresholding. The returned relaxed.projector is the primary convex estimate; loadings are its leading rank eigenvectors and are not claimed to be individually identified across repeated eigenvalues.

Usage

fantope_pca(
  x,
  rank,
  tau,
  center = c("mean", "none"),
  scale = FALSE,
  covariance_divisor = c("n", "n-1"),
  rho = 1,
  solver_tol = 1e-07,
  solver_max_iter = 5000L,
  symmetry_tol = sqrt(.Machine$double.eps),
  strict = TRUE,
  keep_operator = TRUE
)

Arguments

Value

A fantope_pca_fit inheriting from hd_pca_fit.

References

Vu, V. Q., Cho, J., Lei, J., and Rohe, K. (2013). Fantope projection and selection: a near-optimal convex relaxation of sparse PCA. NeurIPS 26.

Examples

x <- rbind(c(3, 0), c(-3, 0), c(0, 1), c(0, -1))
fantope_pca(x, rank = 1, tau = 0.1)

Feng-Jiang-Liu-Xiong max-sum panel-independence test

Description

Implements the three procedures in Feng, Jiang, Liu and Xiong (2022) for serially uncorrelated panel errors. For OLS residual correlations \hat\rho_{ij},

S_N=\sum_{i<j}T\hat\rho_{ij}^2,\qquad L_N=\max_{i<j}|\hat\rho_{ij}|.

The sum component is (S_N-\mu_N)/N, where

\mu_N=\frac{T}{(T-p)^2}\sum_{i<j}\mathrm{tr}(P_iP_j),

and the max component is TL_N^2-4\log N+\log\log N. Their primary max-sum statistic is C_N=\min(p_L,p_S) with calibrated p-value 2C_N-C_N^2.

Usage

feng_jiang_liu_xiong_panel_independence_test(
  panel,
  regressors = NULL,
  component = c("max-sum", "max", "sum"),
  keep_correlations = FALSE
)

Arguments

Value

An object inheriting from htest, with all three component statistics and p-values retained in components.

References

Feng, L., Jiang, T., Liu, B. and Xiong, W. (2022). Max-Sum Tests for Cross-Sectional Independence of High-Dimensional Panel Data. Annals of Statistics, 50, 1124-1143. doi:10.1214/21-AOS2142

Examples

panel <- matrix(c(-2, 1, 0, 2, -1, 3, 1, -2,
                  2, 1, -3, 1, 3, -1, 2, -2), 4, 4)
feng_jiang_liu_xiong_panel_independence_test(panel)

Feng–Lan–Liu–Ma high-dimensional maximum alpha test

Description

Computes the maximum squared unrestricted OLS intercept t statistic, M=\max_i t_i^2, and its centered value M-2\log N+\log\log N. The limiting cdf is \exp\{-\pi^{-1/2}\exp(-x/2)\}. The upper tail is evaluated without subtracting a cdf numerically.

Usage

feng_lan_liu_ma_alpha_max_test(returns, factors = NULL)

Arguments

Value

An asymptotic Gumbel-calibrated upper-tail alpha-test object.

References

Feng, L., Lan, W., Liu, B., and Ma, Y. (2022). High-dimensional test for alpha in linear factor pricing models. Journal of Econometrics. doi:10.1016/j.jeconom.2021.07.011.

Examples

f <- cbind(seq(-1, 1, length.out = 12))
y <- outer(seq_len(12), 1:3, function(i, j) cos(i + 2 * j))
feng_lan_liu_ma_alpha_max_test(y, f)

Feng–Liu–Ma high-dimensional white-noise test

Description

Implements the feasible max statistic, the diagonal-deleted U-statistic sum test, and their Fisher combination from Feng, Liu and Ma. The primary paper uses n (not n - h) in every lagged sample covariance and uses the ordered-pair denominator n * (n - 1) in both the sum statistic and the feasible trace estimate.

Usage

feng_liu_ma_white_noise_test(
  x,
  lag = 1L,
  component = c("fisher", "sum", "max"),
  center = c("none", "mean"),
  keep_lag = FALSE
)

Arguments

Value

An htest object with raw feasible components and diagnostics.

References

Feng, L., Liu, B. and Ma, Y. Testing for High-Dimensional White Noise. Statistica Sinica. doi:10.5705/ss.202023.0300.

Examples

t <- seq_len(18)
x <- cbind(sin(t), cos(t / 2), sin(t / 3 + 0.2))
feng_liu_ma_white_noise_test(x, lag = 2)

Feng–Liu spatial-rank sphericity tests

Description

Implements the Spearman- and Kendall-type high-dimensional sphericity tests of Feng and Liu (2017). method = "spearman" uses

\widehat{\operatorname{tr}(\Omega^2)}= \{2n(n-1)(n-2)(n-3)\}^{-1} \sum^*(U_{ij}^\mathsf{T}U_{kl})(U_{kj}^\mathsf{T}U_{il})

and \widetilde Q=4p\widehat{\operatorname{tr}(\Omega^2)}-1. method = "kendall" removes the factor two, replaces the summand by (U_{ij}^\mathsf{T}U_{kl})^2, and uses \widetilde Q=p\widehat{\operatorname{tr}(\Omega^2)}-1. The star means that all four ordered indices are distinct. Both divide by the same \sigma_0 used by the spatial-sign test.

Usage

feng_liu_rank_sphericity_test(
  x,
  method = c("spearman", "kendall"),
  alpha = 0.05,
  keep_pair_signs = FALSE
)

Arguments

Details

Pair directions are computed once in C++. For each unordered four-set, its 24 ordered terms are reduced exactly to the three disjoint pairings. A tied pair has the literal direction U(0) = 0, is counted in diagnostics, and is never dropped from a denominator.

Value

An object of class c("hd_sphericity_test", "htest") containing the exact ordered sum, denominator, trace estimate, pair-tie diagnostics, and optional pair signs.

References

Feng, L. and Liu, B. (2017). High-dimensional rank tests for sphericity. Journal of Multivariate Analysis, 155, 217–233. doi:10.1016/j.jmva.2017.01.005.

Examples

x <- rbind(c(-2, 0), c(-1, 1), c(0, -2), c(1, 2), c(3, -1), c(2, 1))
feng_liu_rank_sphericity_test(x, method = "kendall")

Feng–Zhang–Liu spatial-rank proportionality test

Description

Implements the high-dimensional spatial-rank test of Feng, Zhang, and Liu (2022). For each group, the within component averages \{U(X_i-X_j)^T U(X_k-X_l)\}^2 over four ordered, mutually distinct indices. The cross component averages over an ordered unequal pair in each group. If the corresponding unscaled averages are A_1,A_2,C_{12}, then

T_{HT}=p(A_1+A_2-2C_{12}).

Usage

feng_spatial_rank_proportionality_test(
  x,
  y,
  alpha = 0.05,
  zero_tol = 0,
  strict = TRUE
)

Arguments

Details

The feasible null variance in the primary paper is implemented literally:

\widehat\sigma_{0,n}^2= \frac{4(n_1^{-1}+n_2^{-1})^2}{(n_1+n_2)(p+2)^2} p^2(n_1A_1+n_2A_2).

The statistic is calibrated by the upper normal tail. Pairwise ties use the paper's convention U(0)=0. A non-positive feasible variance is not floored or replaced.

The primary paper normalizes its shape matrix \Lambda to trace one, and the estimand is p\,tr\{(S_1-S_2)^2\} for the SSCM/Kendall functionals. This is not the trace-p shape-scale formula printed in the current book draft; multiplying \Lambda by p requires the corresponding powers of p in every variance trace.

Value

An hd_proportionality_test object. Both raw U-statistic averages and their p-scaled contributions are returned.

References

Feng, L., Zhang, X., and Liu, B. (2022). High-dimensional proportionality test of two covariance matrices and its application to gene expression data. Statistical Theory and Related Fields, 6, 161–174. doi:10.1080/24754269.2021.1984373.

Examples

x <- matrix(seq_len(40), 8, 5) + matrix(c(-1, 0, 1, 0), 8, 5)
y <- matrix(seq_len(45), 9, 5) + matrix(c(0, 1, -1), 9, 5)
feng_spatial_rank_proportionality_test(x, y)

Feng–Sun scalar-invariant one-sample spatial-sign test

Description

Tests H_0:\theta=\mu_0 against an unrestricted location alternative with the scalar-invariant high-dimensional spatial-sign procedure of Feng and Sun (2016). For every unordered pair i<j, the function jointly fits a leave-two-out location \widehat\theta_{ij} and positive diagonal scale \widehat D_{ij} from all observations except i,j. The raw statistic is

T_{SS}=\frac{2}{n(n-1)}\sum_{i<j} U\{\widehat D_{ij}^{-1/2}(X_i-\mu_0)\}^\mathsf{T} U\{\widehat D_{ij}^{-1/2}(X_j-\mu_0)\}.

In particular, \widehat\theta_{ij} does not enter this numerator.

Usage

feng_sun_one_sample_test(x, mu = NULL, tol = 1e-07, max_iter = 500L)

Arguments

Details

The implementation also supplies the feasible null calibration from the original paper, which is needed for an operational test but is omitted from the short presentation in the book chapter:

\widehat{\operatorname{tr}(R^2)}= \frac{p^2}{n(n-1)}\sum_{i\ne j} \left[U\{\widehat D_{ij}^{-1/2}(X_i-\widehat\theta_{ij})\}^\mathsf{T} U\{\widehat D_{ij}^{-1/2}(X_j-\widehat\theta_{ij})\}\right]^2,

\widehat\sigma_n^2= \frac{2\widehat{\operatorname{tr}(R^2)}}{n(n-1)p^2}, \qquad Z_{SS}=T_{SS}/\widehat\sigma_n.

The p-value is the upper tail of the asymptotic standard normal law. The vector alternative is conventionally labelled two.sided, while large positive quadratic evidence is the rejection direction.

Each leave-out fit follows the paper's diagonal HR recursion. Its initial location and scale are the leave-out sample mean and marginal sample variances. Iteration stops only when both the infinity norm of the mean sign equation and the maximum diagonal-scale equation residual are at most tol. Failure to converge, a coincident training observation (for which the inverse-radius location update is undefined), a non-positive marginal scale, or a non-positive feasible variance raises an error. No ridge, absolute-value repair, or numerical floor is used. The paper itself notes that general existence, uniqueness, and convergence of this diagonal HR recursion are not established.

Before fitting, null-centred data are divided by a safe scale separately in every variable. This is an algebraically neutral coordinatewise change of units for the Feng–Sun statistic and protects both very large and very small finite inputs. Returned leaveout.scale.diagonal values are mapped back to input-coordinate geometry and canonically normalised so their largest diagonal entry is one; the scale equations identify the diagonal only up to a common multiplier. The corresponding canonical log diagonals are also returned, retaining information if a display-scale entry underflows. leaveout.location is in the original input units.

The literal algorithm requires n(n-1)/2 separate iterative fits and is therefore computationally intensive. Its compiled implementation is intended for faithful, moderate-sample use rather than silently replacing the published statistic by a full-sample approximation.

Value

An object of class c("hd_location_test", "htest"). Components include T.SS, the ordered-pair trace estimate, sigma2.hat, pairwise numerator and variance contributions, and every leave-two-out location and canonical diagonal-scale fit. Diagnostics report equation residuals, iteration counts, exact zero signs, internal scaling, and the no-repair policy.

References

Feng, L. and Sun, F. (2016). Spatial-sign based high-dimensional location test. Electronic Journal of Statistics, 10, 2420–2434. doi:10.1214/16-EJS1176.

Examples

set.seed(2608)
x <- matrix(stats::rt(48, df = 5), 8, 6)
feng_sun_one_sample_test(x, tol = 1e-6)

Feng–Wang pairwise-difference-quantile spatial-sign test

Description

Tests equality of two high-dimensional elliptical location vectors using the pairwise-difference-quantile (PDQ) spatial-sign method of Feng and Wang (2026). Observations are rows and variables are columns. The method allows arbitrary within-group correlation and unequal scatter matrices.

Usage

feng_wang_pdq_two_sample_test(
  x,
  y,
  quantile_prob = 0.25,
  level = 0.05,
  B = 999L,
  seed = NULL,
  keep_bootstrap = FALSE,
  tol = 1e-08,
  max_iter = 1000L,
  strict = TRUE
)

Arguments

Details

For group k and coordinate j, let M_k=n_k(n_k-1)/2 and form the M_k unordered absolute pairwise differences. The PDQ scale is the exact empirical U-quantile

q_{kj}=\inf\{t:M_k^{-1}\sum_{i<l} 1(|X_{kij}-X_{klj}|\leq t)\geq a\},

where a is quantile_prob. Thus the implementation selects sorted difference number \lceil aM_k\rceil; it does not use an interpolated sample quantile. The diagonal standardizer has entries q_{kj}^2.

After computing each full-sample spatial median under its PDQ standardizer, the observed full-sample statistic is the cross-centred quantity

\widehat R=-{1\over n_1n_2}\sum_{i,j} U\{\widehat D_1^{-1/2}(X_{1i}-\widehat\theta_2)\}^{T} U\{\widehat D_2^{-1/2}(X_{2j}-\widehat\theta_1)\}.

This is not replaced by the fitted-sign quadratic expansion. From fitted within-group signs the function constructs the published matrices \widehat K_1,\widehat K_2,\widehat K_3 and subtracts the empirical diagonal term

\widehat b=\sum_{k=1}^2 n_k^{-2}\sum_i \widehat S_{ki}^{T}\widehat K_k\widehat S_{ki}.

The reported statistic is T_{PDQ}=\widehat R-\widehat b; large values reject the two-sided location null.

Calibration fixes all nuisance estimates. For each bootstrap draw, independent Rademacher multipliers are applied to the fitted within-group signs, the published quadratic form Q^* is evaluated, and T^*=Q^*-\widehat b is recorded. The paper's primary rule rejects strictly when T_{PDQ} exceeds the lower empirical (1-\mathrm{level}) quantile of the bootstrap values. Because the paper does not prescribe a finite-B p-value, the htest p.value uses the explicitly documented conservative convention

p_{MC}={1+\#\{T^*\geq T_{PDQ}\}\over B+1}.

Ties enter the upper tail. The primary strict-quantile decision and the auxiliary p-value decision are both returned and can differ on the finite bootstrap grid.

The spatial median is computed with a modified Weiszfeld recursion and is accepted when its normalized subgradient residual is at most tol. strict = TRUE makes non-convergence an error; strict = FALSE warns and continues from the last iterate when every published downstream quantity remains defined. A fitted zero residual is always an error because the definition of \widehat G_k requires its inverse radius. Zero cross-centred residuals use the paper's convention U(0)=0 and are counted. Zero PDQ scales, singular \widehat G_k, and non-positive diagonal-deleted bootstrap variance fail explicitly. No ridge, scale floor, generalized inverse, absolute-value correction, or random perturbation is used.

A common coordinatewise affine preconditioning is used internally to protect calculations with extreme but finite units. It is algebraically neutral under the method's common translation and common nonzero coordinatewise scaling invariances. Actual input-unit PDQ diagonals are returned when representable; logarithms and a canonical maximum-one version retain scale ratios when squaring would overflow or underflow.

Value

An object of class c("hd_location_test", "htest"). Its components retain the PDQ quantiles and exact order-statistic ranks, both spatial-median fits, fitted and cross-centred signs, the \widehat\Omega_k, \widehat G_k, bridge and K matrices, observed and fitted quadratic decompositions, diagonal-deletion kernels, bootstrap variance, tail count, and optional draws. diagnostics records both finite-B decisions, convergence, zeros, conditioning, numerical representation, seed, and the no-repair policy.

References

Feng, L. and Wang, H. (2026). High-dimensional two-sample test for elliptical symmetry distribution. arXiv:2605.03265. doi:10.48550/arXiv.2605.03265.

Examples

v <- rbind(c(1, 2, 3), c(2, -3, 1), c(-4, 1, 2), c(3, 4, -2))
x <- rbind(v, -v)
y <- rbind(1.1 * v, -1.1 * v) + rep(c(0.2, -0.1, 0.15), each = 8)
feng_wang_pdq_two_sample_test(x, y, B = 19, seed = 2605)

Feng–Zhang–Liu high-dimensional two-sample spatial-rank test

Description

Tests equality of two high-dimensional elliptical location vectors using the leave-two-out spatial-rank statistic of Feng, Zhang, and Liu (2020). Observations are rows and variables are columns. For distinct i,j in group 1 and s,t in group 2, let \widehat D_{(i,j,s,t)} be the weighted pool of the two diagonal spatial-rank scale fits after deleting i,j and s,t, respectively. The raw statistic is

T_n=\{n_1(n_1-1)n_2(n_2-1)\}^{-1} \sum_{i\ne j}\sum_{s\ne t} U\{\widehat D_{(i,j,s,t)}^{-1/2}(X_{1i}-X_{2s})\}^{\mathsf T} U\{\widehat D_{(i,j,s,t)}^{-1/2}(X_{1j}-X_{2t})\}.

Here U(v)=v/\lVert v\rVert for nonzero v, and U(0)=0.

Usage

feng_zhang_liu_spatial_rank_test(
  x,
  y,
  tol = 1e-07,
  max_iter = 500L,
  scale_identification = c("geometric", "paper_trace")
)

Arguments

Details

The p-value uses the feasible null calibration in the original article. Its two within-group trace estimates use all ordered quadruples of mutually distinct observations and their leave-four-out diagonal scale fits. The cross trace uses all ordered within-group pairs, pooled leave-two-out scales, and the published denominator n_1^2n_2^2. If these estimates are denoted by \widehat{\operatorname{tr}(R_1^2)}, \widehat{\operatorname{tr}(R_2^2)}, and \widehat{\operatorname{tr}(R_1R_2)}, then

\widehat\sigma_n^2= \frac{\widehat{\operatorname{tr}(R_1^2)}} {2n_1(n_1-1)p^2}+ \frac{\widehat{\operatorname{tr}(R_2^2)}} {2n_2(n_2-1)p^2}+ \frac{\widehat{\operatorname{tr}(R_1R_2)}}{n_1n_2p^2}.

The p^2 in the second term restores an evident typographical omission in the author TeX: it is required by the symmetric oracle variance and by consistency of the feasible estimator. The reported statistic is Z=T_n/\widehat\sigma_n; large positive values reject, so the p-value is the upper standard-normal tail. The local-alternative oracle variance is not used.

Each diagonal fit follows the published spatial-rank fixed-point update. scale_identification = "paper_trace" applies the article's literal \operatorname{tr}(D)=p normalization. Because separately fitted trace-normalized matrices acquire different scalar representatives after coordinatewise rescaling, their weighted pool is not exactly coordinatewise-scale equivariant in finite samples. The default "geometric" instead imposes unit geometric mean. All leave-out fits then acquire the same harmless scalar under a common coordinatewise rescaling, making the complete statistic exactly coordinatewise-scale invariant while leaving the published fixed-point equation unchanged. Both choices and their distinction are recorded in diagnostics.

At least six observations per group are necessary: a leave-four-out fit must retain at least two observations. More observations may be needed for a particular data set if a retained subset has zero marginal variance or zero spatial-rank energy. Such degeneracy, nonconvergence, and a non-positive feasible variance are errors. No ridge, absolute value, numerical floor, or random perturbation is applied. Pairwise ties follow U(0)=0 and are counted.

The formal asymptotic calibration in the paper assumes a common scatter matrix and its stated high-dimensional trace conditions. The article's numerical study considered unequal scatter, but it did not establish the null limit in that setting; this function therefore does not claim an unequal-scatter guarantee and does not switch to a simulation or bootstrap calibration.

Value

An object of class c("hd_location_test", "htest"). components contains the raw statistic, all three feasible trace and variance terms, exact ordered-sum denominators, block/permutation contributions, and every full/leave-two/leave-four scale fit. diagnostics records convergence, zero directions, preprocessing, identification, the corrected author-TeX factor, applicability, and the no-repair policy.

References

Feng, L., Zhang, X., and Liu, B. (2020). A high-dimensional spatial rank test for the two-sample location problem. Computational Statistics & Data Analysis, 144, 106889. doi:10.1016/j.csda.2019.106889.

Examples

x <- matrix(c(
  0.7, -1.1, 0.2, -0.4, 0.8, 1.2, 1.1, 0.3, -0.9,
  -1.2, -0.5, 0.7, 0.2, 1.4, -0.4, 1.5, -0.2, 0.5
), ncol = 3, byrow = TRUE)
y <- matrix(c(
  -0.3, 0.7, -0.1, 1.2, -0.6, 0.8, -1.1, 0.2, 1.3,
  0.5, 1.1, -0.7, 0.9, -1.3, 0.4, -0.6, -0.4, -1.2
), ncol = 3, byrow = TRUE)
feng_zhang_liu_spatial_rank_test(x, y, tol = 1e-6)

Feng–Zou–Wang two-sample multivariate-sign test

Description

Tests equality of two multivariate elliptical location vectors with the scalar-invariant spatial-sign procedure of Feng, Zou, and Wang (2016). Observations are rows and variables are columns. If (\widehat\theta_{k,i},\widehat D_{k,i}) is the diagonal HR fit in group k after deleting observation i, the raw statistic is

R_n=-\frac{1}{n_1n_2}\sum_{i=1}^{n_1}\sum_{j=1}^{n_2} U\{\widehat D_{1,i}^{-1/2}(X_{1i}-\widehat\theta_{2,j})\}^{\mathsf T} U\{\widehat D_{2,j}^{-1/2}(X_{2j}-\widehat\theta_{1,i})\}.

Thus, each sign uses its own group's leave-one-out scale but the other group's leave-one-out location. This crossed construction removes the high-dimensional location-estimation bias without estimating and subtracting a separate bias term.

Usage

feng_zou_wang_two_sample_sign_test(x, y, tol = 1e-07, max_iter = 500L)

Arguments

Details

The p-value uses the feasible null calibration in Proposition 2 of the original article, not the oracle variance for local alternatives. Let \widetilde U_{ki}=U\{\widehat D_{k,i}^{-1/2} (X_{ki}-\widehat\theta_{k,i})\} and \widehat c_k=n_k^{-1}\sum_i \|\widehat D_{k,i}^{-1/2}(X_{ki}-\widehat\theta_{k,i})\|^{-1}. With \widehat D_1,\widehat D_2 denoting the two full-sample diagonal HR fits, the three trace estimators are

\widehat{\operatorname{tr}(A_1^2)}= \frac{p^2\widehat c_2^2\widehat c_1^{-2}}{n_1(n_1-1)} \sum_{k=1}^{n_1}\sum_{\ell\ne k} (\widetilde U_{1\ell}^{\mathsf T}\widehat D_2^{-1/2} \widehat D_1^{1/2}\widetilde U_{1k})^2,

with the group indices reversed for \widehat{\operatorname{tr}(A_2^2)}, and

\widehat{\operatorname{tr}(A_3^{\mathsf T}A_3)}= \frac{p^2}{n_1n_2}\sum_{\ell=1}^{n_1}\sum_{k=1}^{n_2} (\widetilde U_{1\ell}^{\mathsf T}\widetilde U_{2k})^2.

These are combined as

\widehat\sigma_n^2= \frac{2\widehat{\operatorname{tr}(A_1^2)}}{n_1(n_1-1)p^2}+ \frac{2\widehat{\operatorname{tr}(A_2^2)}}{n_2(n_2-1)p^2}+ \frac{4\widehat{\operatorname{tr}(A_3^{\mathsf T}A_3)}}{n_1n_2p^2}.

The reported statistic is Z=R_n/\widehat\sigma_n; large positive values reject, so the p-value is the upper standard-normal tail.

Every full-sample and leave-one-out fit follows the published diagonal HR recursion, initialized by the applicable sample mean and marginal sample variances. Convergence requires both estimating-equation residuals to be at most tol. A coincident training observation makes the inverse-radius update undefined and is an error. The spatial-sign convention is U(0)=0 in the numerator, but a zero own-sample leave-one-out radius is also an error because its reciprocal is required by \widehat c_k. Nonconvergence, non-positive scales, and non-positive feasible variance are reported without ridge, absolute-value, floor, or perturbation repairs.

Internally, both groups receive one common, safe coordinatewise affine transformation. This is algebraically neutral under the method's shift and coordinatewise nonzero-scaling invariance and protects extreme finite units. Returned input-coordinate diagonal scales are canonically normalized to have largest diagonal entry one; log diagonals retain ratios that underflow in the display-scale version.

The paper recommends a bootstrap calibration when the dimension is small (it gives p\le 50 as a guide) or grows at order n^2 or faster. This function implements the article's feasible asymptotic-normal test; diagnostics flag those finite-sample regimes. It does not silently switch calibration or reproduce the paper's simulation procedure.

Value

An object of class c("hd_location_test", "htest"). In addition to the test result, components contains both full-sample fits, all leave-one-out fits and directions, inverse radii, bridge diagonals, ordered/cross pair contributions, three trace estimates, and all three variance contributions. diagnostics contains iteration counts, equation residuals, minimum training radii, zero-sign counts, numerical scaling information, asymptotic-regime flags, and the no-repair policy.

References

Feng, L., Zou, C., and Wang, Z. (2016). Multivariate-sign-based high-dimensional tests for the two-sample location problem. Journal of the American Statistical Association, 111(514), 721–735. doi:10.1080/01621459.2015.1035380.

Examples

set.seed(2016)
x <- matrix(stats::rt(48, df = 5), 8, 6)
y <- matrix(stats::rt(54, df = 5), 9, 6)
feng_zou_wang_two_sample_sign_test(x, y, tol = 1e-6)

Feng–Zou–Wang–Zhu scale-invariant Behrens–Fisher test

Description

Tests equality of two high-dimensional mean vectors without requiring equal covariance matrices. Observations are rows and variables are columns. Let \gamma=n_1/n_2, let \widehat\sigma_{sk}^2 be the unbiased sample variance of variable k in group s, and define

A_k=(\bar X_{1k}-\bar X_{2k})^2- \widehat\sigma_{1k}^2/n_1-\widehat\sigma_{2k}^2/n_2.

The initial statistic is

T_{BF}=\sum_k \frac{A_k}{\widehat\sigma_{1k}^2+ \gamma\widehat\sigma_{2k}^2}.

Usage

feng_zou_wang_zhu_two_sample_test(x, y)

Arguments

Details

The reported statistic subtracts the Feng–Zou–Wang–Zhu plug-in estimate of the asymptotic null centering and divides by their ratio-consistent null standard deviation. If \widehat\kappa_{sk}=n_s^{-1}\sum_i (X_{sik}-\bar X_{sk})^3 and D_k=\widehat\sigma_{1k}^2+ \gamma\widehat\sigma_{2k}^2, the centering is

\widehat\mu_n=\widehat b_1+\widehat b_2,

with

\widehat b_1=\sum_k\left\{ \frac{2\widehat\sigma_{1k}^4}{n_1(n_1-1)D_k^2}+ \frac{2\gamma\widehat\sigma_{2k}^4}{n_2(n_2-1)D_k^2} \right\},

and

\widehat b_2=\sum_k\frac{2}{D_k^3} \left(\frac{\widehat\kappa_{1k}}{n_1}- \frac{\gamma\widehat\kappa_{2k}}{n_2}\right)^2.

In particular, the second term of \widehat b_1 contains one power of \gamma. This follows the published paper; the current book draft at ch2_location.tex:685 has an extra power of \gamma. The displayed population expansion in that draft is also an asymptotic centering up to a smaller-order remainder, not an exact finite-sample expectation.

The two within-group covariance trace estimates use exactly the paper's leave-four-out marginal variances and denominator 2P_{n_s}^4. The cross trace removes two observations from each group and uses denominator 4P_{n_1}^2P_{n_2}^2. At least six observations are consequently required in each group. The C++ kernel sums the 24 orderings associated with each unordered within-group quadruple from a local 4\times4 dual Gram matrix. This reduces computation without changing the published estimator and never constructs a p\times p matrix.

Both samples are internally translated by a common anchor and divided by a common scale within each variable. This is algebraically neutral and protects the statistic under very large or small measurement units. Every full-sample and leave-out combined marginal variance, and the final null variance estimate, must have the sign required by the original formulas. No ridge, absolute-value repair, or variance floor is applied.

Value

An object of class c("hd_location_test", "htest"). The upper-tail normal p-value is based on Z=(T_{BF}-\widehat\mu_n)/\widehat\sigma_n. components contains T.BF, Q3 = n1 * T.BF, the two centering terms, all three published leave-out trace estimates, their sample-size coefficients and permutation counts, and coordinate-level quantities. diagnostics records the exact leave-out policy and the absence of numerical repair.

References

Feng, L., Zou, C., Wang, Z., and Zhu, L. (2015). Two-sample Behrens–Fisher problem for high-dimensional data. Statistica Sinica, 25, 1297–1312. doi:10.5705/ss.2014.048.

Examples

set.seed(2411)
x <- matrix(rnorm(56), 8, 7)
y <- matrix(rnorm(63, 0.15), 9, 7)
feng_zou_wang_zhu_two_sample_test(x, y)

Fisher–Sun–Gallagher fourth-to-second trace sphericity test

Description

For N=n+1 observations, this implements the primary paper's unbiased estimators \hat a_2 and \hat a_4 of p^{-1}\operatorname{tr}(\Sigma^2) and p^{-1}\operatorname{tr}(\Sigma^4). With c=p/n, the statistic

\sqrt{\frac{np}{8(8+12c+c^2)}} \left(\frac{\hat a_4}{\hat a_2^2}-1\right)

has an upper standard-normal null calibration. The calculation is performed after a common scale normalization, which is algebraically neutral.

Usage

fisher_sun_gallagher_sphericity_test(x)

Arguments

Value

An hd_covariance_test object retaining every trace and coefficient.

References

Fisher, T. J., Sun, X. and Gallagher, C. M. (2010). Journal of Multivariate Analysis, 101, 2554–2570. doi:10.1016/j.jmva.2010.07.004.

Examples

set.seed(43)
x <- matrix(rnorm(60), nrow = 15, ncol = 4)
fisher_sun_gallagher_sphericity_test(x)

Feng–Lan–Liu–Ma adaptive Gaussian alpha test

Description

Combines the Pesaran–Yamagata sum p-value and the Feng–Lan–Liu–Ma max p-value with the primary paper's two-test Bonferroni rule

p_{\rm COM}=\min\{1,2\min(p_{\rm PY},p_{\rm MAX})\}.

It deliberately does not implement the generic Cauchy benchmark introduced in the book draft, which is not the combination proposed in the cited Feng–Lan–Liu–Ma method.

Usage

gaussian_alpha_combination_test(returns, factors = NULL, p0 = 0.1, delta = 1)

Arguments

Value

An alpha-test object containing both complete component tests.

References

Feng, L., Lan, W., Liu, B., and Ma, Y. (2022). High-dimensional test for alpha in linear factor pricing models. Journal of Econometrics. doi:10.1016/j.jeconom.2021.07.011.

Examples

f <- cbind(seq(-1, 1, length.out = 12))
y <- outer(seq_len(12), 1:3, function(i, j) sin(i * j))
gaussian_alpha_combination_test(y, f)

Classical Gaussian likelihood-ratio test for a covariance matrix

Description

Tests H_0:\Sigma=\Sigma_0 with the Gaussian likelihood-ratio statistic

n\{\operatorname{tr}(\Sigma_0^{-1}S_n) -\log|\Sigma_0^{-1}S_n|-p\},

where S_n=n^{-1}\sum_i(X_i-\bar X)(X_i-\bar X)' when center = TRUE. Its fixed-dimensional reference distribution is chi-squared with p(p+1)/2 degrees of freedom. The sample covariance must be positive definite; no pseudo-determinant or ridge is substituted.

Usage

gaussian_covariance_lrt(x, sigma0, center = TRUE)

Arguments

Value

An object of class c("hd_covariance_test", "htest").

References

Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis, 3rd ed. Wiley.

Examples

x <- matrix(rnorm(80), 20, 4)
gaussian_covariance_lrt(x, diag(4))

Gaussian graphical lasso with an off-diagonal penalty

Description

With the centered empirical covariance S_n, solves the Yuan–Lin off-diagonal graphical-lasso program

\min_{\Omega\succ0} \operatorname{tr}(S_n\Omega)-\log\det(\Omega) +\lambda\sum_{i\ne j}|\omega_{ij}|.

Diagonal entries are never penalized. The optimizer uses an SPD-preserving proximal-gradient step with explicit majorization backtracking and returns the full diagonal/off-diagonal subgradient KKT residual. With lambda = 0, the exact inverse is returned only when S_n is strictly positive definite; no pseudoinverse or ridge is substituted.

Usage

gaussian_graphical_lasso(
  x,
  lambda,
  center = TRUE,
  divisor = c("n", "n-1"),
  solver_tol = 1e-07,
  solver_max_iter = 100000L,
  initial_step = 1,
  max_backtracking = 100L,
  strict = TRUE
)

Arguments

Value

A gaussian_precision_fit list. estimate is non-NULL only when SPD, objective-descent, relative-update, and KKT certificates all pass.

References

Yuan, M. and Lin, Y. (2007). Model selection and estimation in the Gaussian graphical model. Biometrika, 94, 19–35. doi:10.1093/biomet/asm018.

Examples

x <- rbind(c(-2, 0), c(-1, -1), c(1, 1), c(2, 0))
gaussian_graphical_lasso(x, lambda = 0.2)

Oracle Gaussian linear discriminant classifier

Description

Constructs the exact Gaussian log-likelihood-ratio score

(z-(\mu_1+\mu_2)/2)'\Omega(\mu_1-\mu_2) +\log(\pi_1/\pi_2).

Exactly one of covariance and precision must be supplied. No ridge, generalized inverse, or positive-definite repair is used.

Usage

gaussian_lda_oracle(
  location1,
  location2,
  covariance = NULL,
  precision = NULL,
  prior = "equal",
  levels = c("class1", "class2"),
  feature_names = NULL
)

Arguments

Value

A valid hd_classifier_fit.

References

Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis, 3rd ed. Wiley.

Examples

fit <- gaussian_lda_oracle(c(1, 0), c(-1, 0), covariance = diag(2))
predict(fit, rbind(c(2, 0), c(-2, 0)))

Full-covariance Gaussian-mixture EM

Description

Implement the Chapter 7 component-specific full-covariance Gaussian EM updates. Means, covariance arrays, and mixing proportions are all required: this function has no hidden start. The E-step uses Cholesky solves and log-sum-exp; no inverse or pseudo-inverse is formed. The M-step covariance has maximum-likelihood divisor sum(tau[, k]).

Usage

gaussian_mixture_em(
  x,
  means,
  covariances,
  proportions,
  ridge = 0,
  solver_tol = 1e-08,
  monotone_tol = 1e-10,
  rank_tol = 1e-10,
  symmetry_tol = 1e-12,
  solver_max_iter = 200L,
  strict = TRUE,
  keep_responsibilities = TRUE
)

Arguments

Details

The strict default ridge = 0 implements ordinary EM. A positive explicit ridge adds ridge * I after every raw covariance M-step and is labelled regularized EM; unpenalized likelihood monotonicity is then not claimed. Singular initial or updated covariance matrices error instead of receiving an eigenvalue floor. Gaussian-mixture likelihood is unbounded without further restrictions, so a numerical EM fixed point is not a global-MLE certificate.

Value

A gaussian_mixture_em_fit object containing fitted parameters, hard labels (posterior ties use the smallest component), likelihood trace, effective masses, raw covariance ranks, and solver certificates.

References

Feng, L. (2026). High-Dimensional Data Analysis for Elliptical Symmetric Distributions, Chapter 7 (book manuscript). The implementation is the explicitly documented standard Gaussian-mixture E/M construction.

Examples

x <- matrix(c(-2.2, -2, -1.8, 1.8, 2, 2.2), ncol = 1)
gaussian_mixture_em(
  x, means = matrix(c(-2, 2), ncol = 1),
  covariances = array(c(0.2, 0.2), c(1, 1, 2)),
  proportions = c(0.5, 0.5), solver_tol = 1e-6
)

Oracle Gaussian quadratic discriminant classifier

Description

Uses the canonical Gaussian log-likelihood ratio

-\tfrac12\log|\Sigma_1|+\tfrac12\log|\Sigma_2| -\tfrac12d_1^2+\tfrac12d_2^2+\log(\pi_1/\pi_2).

The direct quadratic formula displayed later in the book is twice this score. It has the same boundary but is not the same numerical score.

Usage

gaussian_qda_oracle(
  location1,
  location2,
  covariance1 = NULL,
  covariance2 = NULL,
  precision1 = NULL,
  precision2 = NULL,
  prior = "equal",
  levels = c("class1", "class2"),
  feature_names = NULL
)

Arguments

Value

A valid hd_classifier_fit with a canonical-log-ratio score.

References

Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis, 3rd ed. Wiley.

Examples

fit <- gaussian_qda_oracle(
  c(1, 0), c(-1, 0), covariance1 = diag(2),
  covariance2 = diag(c(2, 1))
)
predict(fit, rbind(c(1, 0), c(-1, 0)), type = "score")

Gaussian Wilks block-independence test

Description

Computes the classical Gaussian likelihood-ratio test for independence of two vector blocks. If \hat\rho_1,\ldots,\hat\rho_m are the sample canonical correlations, m=\min(p,q), then

\Lambda=\prod_{j=1}^m(1-\hat\rho_j^2)

and the Bartlett statistic is

-\{n-1-(p+q+1)/2\}\log\Lambda,

calibrated against \chi^2_{pq}. This is a fixed-dimensional, Gaussian benchmark; it is not a high-dimensional repair of Wilks' test.

Usage

gaussian_wilks_independence_test(x, y)

Arguments

Value

An object inheriting from htest. components contains Wilks' Lambda, its log value, canonical correlations, and the three sample covariance blocks.

References

Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis, 3rd edition, Chapter 8. Wiley.

Examples

x <- matrix(c(-2, 0, -1, 2, 0, -1, 1, 1, 2, -2, 3, 0), 6, 2)
y <- matrix(c(1, -2, 0, 2, -1, 3), 6, 1)
gaussian_wilks_independence_test(x, y)

Generalized spatial-sign principal component analysis

Description

Forms the generalized spatial-sign covariance matrix

n^{-1}\sum_i g(x_i-T)g(x_i-T)',\qquad g(t)=\xi(\lVert t\rVert_2)t,

and decomposes it. Available radial multipliers are "winsor", "quadratic", "ball", "shell", "linear_redescending", "spatial" (ordinary spatial signs), and "identity" (the centered second moment).

Usage

generalized_sign_pca(
  x,
  rank = NULL,
  weight = c("winsor", "quadratic", "ball", "shell", "linear_redescending", "spatial",
    "identity"),
  center = c("kstep_lts", "spatial", "mean", "none"),
  cutoff = c("median_mad", "original_h_order", "user"),
  cutoffs = NULL,
  q1_lower_bound = c("zero", "error"),
  zero_mad = c("limit", "error"),
  lts_steps = 2L,
  divisor = c("n", "nonzero"),
  zero_action = c("zero", "error"),
  tol = 1e-08,
  max_iter = 500L,
  zero_tol = 0,
  keep_operator = TRUE
)

Arguments

Details

The default cutoff convention is the ordinary median/raw-MAD construction stated in Chapter 6 of the book. Because the cited 2024 manuscript was not available for formula verification, "median_mad" is documented here as a book-defined convention and is not claimed to reproduce that manuscript. "original_h_order" implements the 2019 order statistic with h=\lfloor(n+p+1)/2\rfloor and fails when h > n. "user" requires all four named cutoffs. Ball includes Q2; Shell includes both Q1 and Q3. When the MAD is zero, the default exact limiting rule uses coincident cutoffs and never perturbs a denominator. The "spatial" and "identity" families use no cutoffs and reject an explicitly supplied cutoff or cutoffs argument.

For a general elliptical representation X-\mu=RAU, U=Z/\lVert Z\rVert, and S=\sum_l\lambda_l Z_l^2, the population generalized eigenvalue contains the radial variable:

E\{K^2(|R|\sqrt{S}/\lVert Z\rVert) \lambda_j Z_j^2/S\}.

The simpler expression without R is therefore not claimed outside the Gaussian representation. Eigenvector and ordering results require the conditions stated in Chapter 6 for this book-defined construction, not merely a measurable weight.

Value

An object inheriting from hd_pca_fit. eigenvalues contains the full spectrum and loadings contains the requested leading vectors.

References

Raymaekers, J. and Rousseeuw, P. J. (2019). A generalized spatial sign covariance matrix. Journal of Multivariate Analysis, 171, 94–111.

Chapter 6 cites Leyder, S., Raymaekers, J., and Verdonck, T. (2024), Generalized spherical principal component analysis, Statistics and Computing, 34, 104. That manuscript was unavailable for formula verification; the "median_mad" construction above is therefore attributed to the book rather than asserted as a reproduction of the cited article.

Examples

x <- rbind(c(3, 0), c(-3, 0), c(0, 1), c(0, -1))
generalized_sign_pca(
  x, rank = 1, weight = "winsor", center = "none"
)

Generic weighted Hettmansperger–Randles location estimator

Description

Solves the formula-complete weighted location equation in Chapter 2 while retaining the unweighted diagonal Hettmansperger–Randles (HR) scale equation. At iteration m, let

e_i^{(m)}=(D^{(m)})^{-1/2}(X_i-\theta^{(m)}),\quad r_i^{(m)}=\|e_i^{(m)}\|,\quad U_i^{(m)}=e_i^{(m)}/r_i^{(m)}.

The updates are

\theta^{(m+1)}=\theta^{(m)}+(D^{(m)})^{1/2} \frac{\sum_iK(r_i^{(m)})U_i^{(m)}} {\sum_iK(r_i^{(m)})/r_i^{(m)}}

and

D^{(m+1)}=pD^{(m)}\operatorname{diag} \{n^{-1}\sum_iU_i^{(m)}U_i^{(m)T}\}.

Thus weights enter the location equation only; the scale update is always the unweighted HR update printed in the manuscript.

Usage

generic_weighted_hr_location(
  x,
  K = "constant",
  power = 0,
  initial_location = NULL,
  initial_diagonal = NULL,
  tol = 1e-08,
  max_iter = 500L,
  zero_tol = 0,
  keep_history = FALSE
)

Arguments

Details

K may be an R function or one of the auditable built-ins "constant", "inverse_norm", and "power". A user callback is called separately on each scalar radius and must return exactly one finite numeric scalar. The recursion rejects zero radii, a non-positive weighted denominator, non-positive scale iterates, and non-convergence. It never floors a radius, caps a weight, flips a denominator sign, adds a ridge, or returns an unconverged last iterate.

Value

A list of class generic_weighted_hr_location containing the fitted location and diagonal, final directions/radii/weights, equation residuals, and strict convergence diagnostics.

References

Feng, L., Liu, B. and Ma, Y. (2021). An inverse norm sign test for location parameters in high-dimensional data. Journal of Business & Economic Statistics 39, 807–815.

Examples

x <- matrix(c(-2, 1, 0, 3, -1, 2, 1, -3, 2, 0, 4, -2), ncol = 2)
generic_weighted_hr_location(x, K = "constant", tol = 1e-6)

Generalized quadratic discriminant analysis for elliptical populations

Description

Fits the Bose–Pal–Saha Ray–Nayak generalized QDA rule using classical class moments or two explicitly supplied certified location–scatter fits. The rule assigns an observation to class 1 when

\Delta_d^2(x) \ge c\log\{|S_1|/|S_2|\},\qquad 0\le c\le1.

A supplied c is used directly. Otherwise the resubstitution selector (or an explicitly requested leakage-free classical cross-validation selector) enumerates both endpoints, every legal decision breakpoint, and adjacent midpoints. This direct loss calculation remains correct for positive, negative, and zero signed log-determinant contrasts.

Usage

gqda_classifier(
  x,
  y,
  class_fits = NULL,
  c = NULL,
  selection = base::c("resubstitution", "fixed", "cross_validation"),
  folds = 5L,
  seed = NULL,
  strict = TRUE
)

Arguments

Details

This is an equal-prior classifier. The current book's HR-QDA subsection writes an ordinary Gaussian plug-in QDA score, but the cited Yan–Feng–Zhang primary applies the generalized threshold above. No simulation tuning grid, ridge, determinant absolute value, or covariance repair is used.

Value

An hd_classifier_fit. Package scores are positive for class 1.

References

Bose, S., Pal, A., SahaRay, R., and Nayak, J. (2015). Generalized quadratic discriminant analysis. Pattern Recognition, 48(8), 2676–2684. doi:10.1016/j.patcog.2015.02.016.

Examples

x <- rbind(c(-2, 0), c(-1, 1), c(-1, -1), c(-2, 1),
           c(2, 0), c(1, 1), c(1, -1), c(2, -1))
y <- factor(rep(c("left", "right"), each = 4),
            levels = c("left", "right"))
fit <- gqda_classifier(x, y, c = 0)
predict(fit, x)

Gibbons–Ross–Shanken exact alpha test

Description

Tests whether all intercepts are zero in an unconditional linear factor-pricing model. Rows of returns and factors are observations and columns are assets and factors. The joint regression contains both an intercept and all factor columns. With unrestricted OLS residual matrix E, Vhat = E' E / T, and h the residual from regressing the all-ones vector on the factors, the implemented primary statistic is

[(T-N-K)/N][h'h/T]\hat\alpha'\widehat V^{-1}\hat\alpha,

with the exact F distribution having N and T-N-K degrees of freedom under the Gaussian GRS assumptions.

Usage

grs_alpha_test(returns, factors = NULL)

Arguments

Details

This differs from two finite-sample expressions in the book draft. Its displayed slope estimator omits the intercept although its residual formula includes one. It also divides the residual covariance by T-K-1 while retaining the primary statistic scaling for a divisor-T covariance, missing the compensating factor T/(T-K-1). This function follows GRS exactly and reports both divisors for audit.

Value

An object of class c("hd_alpha_test", "htest") with the exact statistic, p-value, fitted intercepts, and OLS diagnostics.

References

Gibbons, M. R., Ross, S. A., and Shanken, J. (1989). A test of the efficiency of a given portfolio. Econometrica, 57, 1121–1152. doi:10.2307/1913625.

Examples

f <- cbind(seq(-1, 1, length.out = 10))
y <- cbind(0.2 + f[, 1] + sin(seq_len(10)),
           -0.1 - 0.5 * f[, 1] + cos(seq_len(10)))
grs_alpha_test(y, f)

Hallin–Paindaveine signed-rank test for elliptical shape

Description

For null shape V_0, the observations are sphericized, their radii are ranked, and their directions are weighted by a score K. The statistic is Hallin and Paindaveine's equation (4.3),

\frac{p(p+2)}{2nE\{K^2(U)\}} \sum_{i,j}K(R_i/(n+1))K(R_j/(n+1)) \{(U_i'U_j)^2-1/p\}.

It has an asymptotic chi-squared reference with p(p+1)/2-1 degrees of freedom. Scores are "sign" (constant), "wilcoxon" (power one), "spearman" (power two), and "vdw" (chi-squared normal scores). Non-sign scores require distinct radii, as in the paper's continuous elliptical model. No random tie breaking is used.

Usage

hallin_paindaveine_shape_test(
  x,
  shape0 = NULL,
  center = NULL,
  score = c("sign", "wilcoxon", "spearman", "vdw"),
  tol = 1e-08,
  max_iter = 500L,
  zero_tol = 0
)

Arguments

Value

An hd_covariance_test object. components retains scores, ranks, signs, radii on the log scale, and the weighted sign scatter matrix.

References

Hallin, M. and Paindaveine, D. (2006). Annals of Statistics, 34, 2707–2756. doi:10.1214/009053606000000731.

Examples

set.seed(38)
x <- matrix(rnorm(30), nrow = 15, ncol = 2)
hallin_paindaveine_shape_test(
  x, diag(2), center = c(0, 0), score = "sign"
)

Construct a classifier from a scoring function

Description

This low-level constructor creates a two-class hd_classifier_fit from a user-supplied scoring function. The function must accept a finite numeric matrix with observations in rows and return one finite score per row. Positive scores select the first class. Package methods use transparent structured linear or quadratic score models instead of this custom route.

Usage

hd_classifier_fit(
  score,
  n_features,
  levels = c("class1", "class2"),
  feature_names = NULL,
  method = "Custom two-class classifier",
  score_scale = "method_threshold",
  tie = c("class1", "class2")
)

Arguments

Value

An object of class hd_classifier_fit.

References

Feng, L. (2026). High-Dimensional Data Analysis for Elliptically Symmetric Distributions, Chapter 5 (book manuscript). This low-level constructor contract is package infrastructure, not a distinct method.

Examples

fit <- hd_classifier_fit(function(z) z[, 1] - z[, 2], 2)
predict(fit, rbind(c(2, 1), c(0, 1)))

High-dimensional Hettmansperger–Randles estimation

Description

Implements Algorithm 2 of Yan, Feng, and Zhang (2025). Starting from the sample spatial median and a supplied positive-definite pilot precision, the method jointly updates location and trace-normalized shape. At iteration k, with e_i^{(k)}=\Sigma_k^{-1/2}(X_i-\mu_k), it applies

\mu_{k+1}=\mu_k+\Sigma_k^{1/2} \frac{n^{-1}\sum_i U(e_i^{(k)})} {n^{-1}\sum_i\|e_i^{(k)}\|^{-1}}

and

\Sigma_{k+1}\ \mathrel{\propto}\ \Sigma_k^{1/2}\mathcal B_h\left\{ n^{-1}\sum_i U(e_i^{(k)})U(e_i^{(k)})^T\right\} \Sigma_k^{1/2},\qquad \operatorname{tr}(\Sigma_{k+1})=p.

Usage

high_dimensional_hr(
  x,
  pilot_precision,
  bandwidth = 3L,
  tol = 1e-08,
  max_iter = 1000L,
  median_tol = 1e-08,
  median_max_iter = 1000L,
  zero_tol = 0,
  scale_estimator = c("paper_qda", "none"),
  strict = TRUE
)

Arguments

Details

The paper sets h=3; the formula permits 0\le h<p. When the default is not explicitly supplied and p\le3, it resolves to min(3, p - 1) and records that value. The pilot's common scalar is not identified by spatial signs, so its inverse is normalized to trace p before iteration.

The primary estimator is a shape estimator, not a covariance-scale estimator. With scale_estimator = "paper_qda", this function additionally implements the paper's QDA scale add-on

\widehat{\operatorname{tr}(\Xi)}= \{\sum_i\|X_i\|^2-n\|\bar X\|^2\}/(n-1),

returning scatter.scale = trace.hat / p and scatter = scatter.scale * shape. This second-moment scale is not part of the robust sign fixed point and requires its finite double representation.

The paper says to repeat until convergence but does not prescribe a norm. Here, convergence means that the maximum of the relative location and Frobenius shape-map updates is at most tol. The spatial-sign score and shape-map residual are returned separately; they are diagnostics rather than an independently claimed estimating-equation certificate.

A supplied spatial_sign_precision() fit is accepted only when its feasibility/KKT certificate made it a valid fit. A raw matrix must be symmetric positive definite. Hard banding itself can destroy positive definiteness. Such a failure, a singular pilot, an exact zero standardized residual, overflow, or nonconvergence is reported without a ridge, eigenvalue floor, pseudoinverse, jitter, or perturbation. With strict = FALSE, every such method failure returns estimate = NULL and valid = FALSE; the last iterate appears only under diagnostics.

The book's fixed-pilot score, post-hoc banded raw shape, second banding of the inverse, and HR-centered divisor-n scale are not Algorithm 2. Likewise, the book's displayed r_n+h^{-\alpha} theorem is a review synthesis rather than a finite-sample guarantee stated in this primary.

Value

An object of class high_dimensional_hr_fit. A valid object returns robust location, the final trace-p shape, the unbanded raw.shape fixed-point map, precision = solve(shape), the final SSCMs, optional covariance scatter.scale and scatter, and detailed median, iteration, score, equation, positive-definiteness, reciprocal-condition, scale-identification, and pilot-certificate diagnostics.

References

Yan, G., Feng, L., and Zhang, X. (2025). High-Dimensional Hettmansperger-Randles Estimator and its Applications. arXiv:2505.01669. https://arxiv.org/abs/2505.01669

Examples

x <- rbind(
  c(2, 0), c(-2, 0), c(0, 1), c(0, -1),
  c(1, 1), c(-1, -1), c(1, -1), c(-1, 1)
)
high_dimensional_hr(x, diag(2), bandwidth = 0)

Classical Hotelling location tests

Description

Performs the exact Gaussian one- or two-sample Hotelling T^2 test. Observations are rows and variables are columns. The two-sample test uses the pooled covariance matrix and therefore assumes a common population covariance matrix.

Usage

hotelling_one_sample_test(x, mu = NULL, tol = sqrt(.Machine$double.eps))

hotelling_two_sample_test(x, y, tol = sqrt(.Machine$double.eps))

Arguments

Details

The covariance system is solved by a Cholesky factorization. A generalized inverse is deliberately not used: singular or numerically ill-conditioned covariance matrices invalidate the exact F calibration and produce an error directing the user to a high-dimensional test.

Value

An object of classes hd_location_test and htest. In addition to the standard htest fields, it contains the raw Hotelling T2, the covariance estimate, exact null-distribution metadata, and numerical diagnostics including the reciprocal condition number and solver.

References

Hotelling, H. (1931). The generalization of Student's ratio. Annals of Mathematical Statistics, 2, 360–378.

Examples

set.seed(11)
x <- matrix(rnorm(60), 20, 3)
hotelling_one_sample_test(x)

y <- matrix(rnorm(75, 0.25), 25, 3)
hotelling_two_sample_test(x, y)

Hettmansperger-Randles affine-equivariant location and shape

Description

Alternates the location and Tyler-type shape updates stated in Chapter 1. The returned shape has trace equal to the data dimension.

Usage

hr_estimator(
  x,
  initial_location = NULL,
  initial_shape = NULL,
  tol = 1e-08,
  max_iter = 500L,
  zero_tol = 0,
  warn = TRUE
)

Arguments

Value

An object of class hd_hr containing location, shape, and convergence diagnostics. In one dimension the method is defined by the conventional sample median and unit shape. This univariate extension does not report the multivariate joint shape-equation residual, which is returned as NA.

References

Hettmansperger, T. P. and Randles, R. H. (2002). A practical affine equivariant multivariate median. Biometrika, 89, 851-860.

Examples

set.seed(5)
x <- matrix(rt(240, df = 4), 60, 4)
hr_estimator(x)

High-dimensional HR generalized quadratic discriminant analysis

Description

Fits the primary Yan–Feng–Zhang classwise high-dimensional HR estimator, including its separate QDA second-moment scale, and plugs the two certified fits into generalized QDA. This is not the ordinary Gaussian HR-QDA score currently displayed in the book.

Usage

hr_gqda(
  x,
  y,
  pilot_precision,
  bandwidth = 3L,
  c = NULL,
  selection = base::c("resubstitution", "fixed", "cross_validation"),
  folds = 5L,
  seed = NULL,
  tol = 1e-08,
  max_iter = 1000L,
  median_tol = 1e-08,
  median_max_iter = 1000L,
  zero_tol = 0,
  strict = TRUE
)

hr_qda(...)

Arguments

Value

An hd_classifier_fit containing both certified HR fits.

References

Yan, G., Feng, L., and Zhang, X. (2025). High-Dimensional Hettmansperger–Randles Estimator and its Applications. arXiv:2505.01669. https://arxiv.org/abs/2505.01669

Examples

x <- rbind(c(-2, 0), c(-1, 1), c(-1, -1), c(-2, 1),
           c(2, 0), c(1, 1), c(1, -1), c(2, -1),
           c(-1.5, .4), c(1.5, -.4), c(-1.7, -.3), c(1.7, .3))
y <- factor(rep(c("left", "right"), each = 6))
hr_gqda(x, y, pilot_precision = diag(2), bandwidth = 0, c = 0)

Influential features PCA clustering

Description

Implements influential features PCA (IF-PCA): standardize every feature with its n - 1 sample standard deviation, rank features by

\psi_{n,j}=\sqrt n\sup_t|\widehat F_{n,j}(t)-\Phi(t)|,

optionally translate and rescale the scores by an empirical null, select features, compute the leading K - 1 left singular vectors, and cluster their rows.

Usage

if_pca(
  x,
  K,
  selection = c("fixed", "hct"),
  threshold = NULL,
  empirical_null = c("mean_sd", "median_mad", "none"),
  ks_convention = c("paper", "software"),
  hct_convention = c("paper", "software"),
  null_scores = NULL,
  null_cdf = NULL,
  null_reps = NULL,
  seed = NULL,
  truncate = FALSE,
  rank_tol = sqrt(.Machine$double.eps),
  kmeans_tol = 1e-10,
  kmeans_max_iter = 100L
)

Arguments

Details

ks_convention = "paper" evaluates the displayed statistic on the n - 1 standardized data. "software" reproduces the distinct official MATLAB score convention by dividing those standardized values once more by sqrt(1 - 1/n), which is equivalent to using divisor n inside the KS calculation. This distinction does not change the PCA matrix. The book's printed n-divisor standardization is not used because the primary method and software both form the PCA matrix with the n - 1 SD.

With fixed selection, threshold is applied to the adjusted scores using a greater-than-or-equal rule. HCT requires either null_scores already on the same final scale, a supplied null_cdf, or explicit null_reps plus seed. A standard KS CDF is not silently substituted because centering and scale are estimated. The original software used 100p null replicates and the paper's numerical study used 2000p; neither computational choice is hidden as a default.

Value

An if_pca_fit object with labels, selected indices and names, standardized data, raw and adjusted KS scores, the embedding, deterministic k-means diagnostics, and full HCT search details when used.

References

Jin, J. and Wang, W. (2016). Influential features PCA for high dimensional clustering. Annals of Statistics, 44, 2323–2359. doi:10.1214/15-AOS1423.

Examples

x <- cbind(
  c(-3, -2.5, -2, -1.5, 1.5, 2, 2.5, 3),
  c(-1, 0, 1, 0, -1, 0, 1, 0),
  c(0, 1, 0, -1, 0, 1, 0, -1)
)
if_pca(x, K = 2, selection = "fixed", threshold = 0,
       empirical_null = "none")

Diagonal-covariance independence classifier

Description

Fits the Gaussian independence rule after replacing the common covariance by its pooled diagonal. variance_divisor = "unbiased" uses n - 2; "mle" uses n. Zero marginal variances are never floored. They either cause a failure or are dropped explicitly.

Usage

independence_classifier(
  x,
  y,
  prior = "equal",
  variance_divisor = c("unbiased", "mle"),
  zero_variance = c("error", "drop"),
  strict = TRUE
)

Arguments

Value

An hd_classifier_fit.

References

Bickel, P. J. and Levina, E. (2004). Some theory for Fisher's linear discriminant function, naive Bayes, and some alternatives when there are many more variables than observations. Bernoulli, 10, 989–1010.

Examples

x <- rbind(c(2, 1), c(1, 2), c(3, 0),
           c(-2, -1), c(-1, -2), c(-3, 0))
independence_classifier(x, rep(c("A", "B"), each = 3))

Feng–Liu–Ma one-sample inverse norm sign test

Description

Tests H_0:\theta=\mu_0 against an unrestricted multivariate location alternative with the inverse norm sign test (INST) of Feng, Liu, and Ma (2021). For every unordered pair i<j, a location and positive diagonal scale are jointly fitted from the sample with observations i,j removed. Writing the fitted scale as \widehat D_{ij}, r_{ij,k}=\lVert\widehat D_{ij}^{-1/2}(X_k-\mu_0)\rVert, and U_{ij,k}=U\{\widehat D_{ij}^{-1/2}(X_k-\mu_0)\}, the primary raw statistic is

T_{INST}=\frac{2}{n(n-1)}\sum_{i<j} r_{ij,i}^{-1}r_{ij,j}^{-1}U_{ij,i}^{\mathsf T}U_{ij,j}.

The pair-specific fitted location is used only to estimate \widehat D_{ij}. It is deliberately not subtracted from either endpoint in this statistic.

Usage

inst_one_sample_test(
  x,
  mu = NULL,
  alpha = 0.05,
  tol = 1e-08,
  max_iter = 500L,
  zero_tol = 0,
  strict = TRUE
)

Arguments

Details

The operational calibration is the direct feasible variance in Section S.3 of the paper's official supplement. If \widetilde\mu_{ij}=(n-2)^{-1}\sum_{k\ne i,j}U_{ij,k}, then

\widehat\sigma_n^2=2n^{-4}\sum_{i\ne j} r_{ij,i}^{-2}r_{ij,j}^{-2} \{(U_{ij,i}-\widetilde\mu_{ij})^{\mathsf T}U_{ij,j}\} \{(U_{ij,j}-\widetilde\mu_{ij})^{\mathsf T}U_{ij,i}\}.

This is an ordered-pair sum with the published 2n^{-4} multiplier; it is not replaced by a finite-sample pair-count denominator. The reported p-value is based only on T_{INST}/\widehat\sigma_n.

For interpretation, the result additionally reports cross-fitted estimates of E(r^{-1}), \nu_{2,IN}=E(r^{-2}), c_{0,IN}=E(r^{-2}), and \operatorname{tr}(R^2). The displayed oracle-factorized variance 2\widehat\nu_{2,IN}^2\widehat{\operatorname{tr}(R^2)}/ \{n(n-1)p^2\} is diagnostic only and never replaces the primary direct variance.

Each diagonal HR recursion starts at the leave-two-out sample mean and marginal sample variances. The common scale from that initialization is retained because the fixed-point equations identify only relative diagonal scales. iteration.stable means that both the relative location update and log-diagonal update are at most tol; the estimating-equation score.residual is returned separately. With strict = TRUE, any required fit that is not stable after max_iter updates is an error. With strict = FALSE, the last iterates are used with a warning and complete diagnostics.

Exact or tolerance-defined zero radii make inverse norm weighting singular. Such a case is an error. A non-finite or non-positive published variance is also an error. No radius perturbation, weight cap, ridge, absolute value, numerical floor, or variance substitution is applied.

Null-centred residuals are internally divided by a positive scale in each coordinate. This is algebraically neutral under a common translation and nonsingular diagonal change of units, and protects the calculation from avoidable overflow. Returned input-coordinate diagonals are canonical displays with largest entry one; their canonical log values retain scale ratio information when a displayed entry underflows. The unnormalised internal diagonals retain the initialization's common scale.

Value

An object of class c("hd_location_test", "htest"). Its components field contains the primary statistic and direct variance, cross-fitted radial moments and trace diagnostics, every pair's kernels, radii, fitted locations, and diagonal scales. Its diagnostics field records iteration stability, score residuals, centring roles, rejection rule, numerical scaling, and the no-repair contract.

References

Feng, L., Liu, B., and Ma, Y. (2021). An inverse norm sign test of location parameter for high-dimensional data. Journal of Business & Economic Statistics, 39, 807–815. doi:10.1080/07350015.2020.1736084.

Feng, L., Liu, B., and Ma, Y. (2020). Supplementary material for An inverse norm sign test of location parameter for high-dimensional data. doi:10.6084/m9.figshare.11914095.v2.

Examples

set.seed(2021)
x <- matrix(stats::rt(40, df = 5), 10, 4)
inst_one_sample_test(x, tol = 1e-6)

Jiang-Wang-Leng direct sparse QDA

Description

Implements the direct quadratic-discriminant estimator of Jiang, Wang and Leng. The interaction estimate minimizes a penalized quadratic trace loss, is symmetrized after optimization, and the linear coefficient minimizes its lasso quadratic loss. This is deliberately not the Dantzig program sometimes misattributed to this paper. The intercept exhaustively checks sorted raw training-score breakpoints and intervening intervals for minimum 0-1 loss.

Usage

jiang_da_qda(
  x,
  y,
  lambda_interaction = NULL,
  lambda_linear = NULL,
  parameter_grid = NULL,
  folds = 5L,
  rho = 1,
  solver_tol = 1e-07,
  solver_max_iter = 10000L,
  strict = TRUE
)

Arguments

Value

An hd_classifier_fit. Strictly positive scores select class 1; exact zero selects class 2, following the primary decision rule.

References

Jiang, B., Wang, X. and Leng, C. (2018). A direct approach for sparse quadratic discriminant analysis. Journal of Machine Learning Research, 19(31), 1-37.

Examples

x <- rbind(
  c(-2, 0), c(-1, 1), c(-1, -1), c(-2, 1),
  c(2, 0), c(1, 2), c(1, -2), c(2, 1)
)
y <- factor(rep(c("left", "right"), each = 4))
fit <- jiang_da_qda(x, y, 0.2, 0.2, solver_tol = 1e-6)

John's classical trace test of Gaussian sphericity

Description

Computes U=p\operatorname{tr}(S^2)/\operatorname{tr}^2(S)-1 and uses mpU/2, with residual degrees of freedom m, as a chi-squared statistic with (p-1)(p+2)/2 degrees of freedom. The multiplier is p; the p+2 multiplier in the book draft belongs to the spatial-sign analogue, not John's covariance statistic.

Usage

john_sphericity_test(x, center = TRUE)

Arguments

Value

An hd_covariance_test object.

References

John, S. (1971). Some optimal multivariate tests. Biometrika, 58, 123–127.

Examples

set.seed(36)
x <- matrix(rnorm(60), nrow = 20, ncol = 3)
john_sphericity_test(x)

K-spatial-median clustering

Description

Alternates exact sample spatial-median center updates with deterministic nearest-center assignments. The fitted objective is the sum of unsquared Euclidean distances. Empty-cluster repair is disabled by default.

Usage

k_spatial_median(
  x,
  K,
  init = "maxmin",
  first_index = 1L,
  tol = 1e-08,
  max_iter = 100L,
  spatial_max_iter = 500L,
  zero_tol = 0,
  empty_action = c("error", "farthest"),
  cycle_action = c("error", "return"),
  ties = "first",
  keep_distances = FALSE
)

Arguments

Value

A k_spatial_median_fit object with labels, centers, objective, and explicit convergence, tie, and repair diagnostics.

References

Zhao, P., Zhuang, D., and Feng, L. (2026). Sparse K-spatial-median clustering for high-dimensional data. arXiv:2605.00598.

Examples

x <- matrix(c(0, 2, 8, 10), ncol = 1)
k_spatial_median(x, K = 2, init = c(1, 4))

Robust factor-number selection from multivariate Kendall eigenvalues

Description

Implements the modified Kendall eigenvalue-ratio (MKER) and transformed contribution-ratio (MKTCR) selectors. With m = min(N, T) for N = p variables and T = n observations, the stabilized eigenvalues are

\widetilde\lambda_j=\lambda_j+c/\sqrt{m}.

MKER maximizes tilde_lambda[j] / tilde_lambda[j + 1]. For MKTCR,

V_j=\sum_{i=j+1}^{m}\widetilde\lambda_i,\qquad \frac{\log(1+\widetilde\lambda_j/V_{j-1})} {\log(1+\widetilde\lambda_{j+1}/V_j)}

is maximized. kmax and strictly positive c are always supplied, and an exact tie selects the smallest index.

Usage

kendall_factor_number(
  x,
  kmax,
  c,
  method = c("mker", "mktcr"),
  zero_tol = 0,
  zero_pair_action = c("error", "zero"),
  eigen_tol = sqrt(.Machine$double.eps)
)

Arguments

Value

A kendall_factor_number object with the selected factor number, raw and stabilized eigenvalues, criterion path, and MKTCR tail sums.

References

Yu, L., He, Y. and Zhang, X. (2019). Robust factor number specification for large-dimensional elliptical factor models. arXiv:1808.09107.

Examples

x <- rbind(c(-3, -1, 0), c(-2, 1, 1), c(-1, -2, 0),
           c(1, 2, 0), c(2, -1, -1), c(3, 1, 0))
kendall_factor_number(x, kmax = 1, c = 0.1, method = "mker")

Multivariate-Kendall principal component analysis

Description

Decomposes the exact multivariate Kendall U-statistic over all unordered pairs. Its default tie convention matches spatial_kendall(): a tied difference contributes the zero matrix while the divisor remains choose(n, 2). The operator is translation invariant; center controls only the reported scores.

Usage

kendall_pca(
  x,
  rank = NULL,
  center = c("mean", "spatial", "none"),
  divisor = c("all_pairs", "nonzero_pairs"),
  ties = c("zero", "error"),
  tol = 1e-08,
  max_iter = 500L,
  zero_tol = 0,
  keep_operator = TRUE
)

Arguments

Value

An object inheriting from hd_pca_fit. eigenvalues contains the full spectrum and loadings contains the requested leading vectors.

References

Han, F. and Liu, H. (2018). ECA: High-dimensional elliptical component analysis in non-Gaussian distributions. Journal of the American Statistical Association, 113, 252–268.

Examples

x <- rbind(c(2, 0), c(-2, 0), c(0, 1), c(0, -1))
kendall_pca(x, rank = 1)

Li–Chen two-sample high-dimensional covariance equality test

Description

Computes the exact location-invariant leave-four-out statistics A_{n_1}, A_{n_2}, and the two-by-two leave-out cross statistic C_{n_1n_2} from Li and Chen (2012). Their combination

T_{LC}=A_{n_1}+A_{n_2}-2C_{n_1n_2}

is unbiased for \operatorname{tr}(\Sigma_1-\Sigma_2)^2. Under the equality null, the primary feasible standard-error estimator is

\hat\sigma_0=2A_{n_1}/n_2+2A_{n_2}/n_1.

The published article labels this quantity once as \hat\sigma_0^2; its units, subsequent ratio-consistency theorem, and rejection rule show that it is the standard error. A non-positive estimate is rejected without an absolute value or floor.

Usage

li_chen_covariance_test(x, y)

Arguments

Value

An hd_covariance_test object retaining original-unit and internally scaled trace estimates.

References

Li, J. and Chen, S. X. (2012). Annals of Statistics, 40, 908–940. doi:10.1214/12-AOS993.

Examples

set.seed(44)
x <- matrix(rnorm(24), nrow = 8, ncol = 3)
y <- matrix(rnorm(27, mean = 0.1), nrow = 9, ncol = 3)
li_chen_covariance_test(x, y)

Li-Shao thresholded sparse QDA

Description

Fits the three-threshold rule of Li and Shao. Class covariances use the maximum-likelihood divisor n_k. Mean differences are retained only when strictly above the first threshold. Entries whose class difference is at most the second threshold are pooled; off-diagonal entries are then kept only when strictly above the third threshold. Diagonals are not removed.

Usage

li_shao_sparse_qda(
  x,
  y,
  threshold_mean = NULL,
  threshold_difference = NULL,
  threshold_covariance = NULL,
  selection = c("specified", "paper_bisection"),
  bisection_tol = NULL,
  ridge = 0,
  strict = TRUE
)

Arguments

Details

Paper bisection uses deterministic leave-one-out endpoint searches. Its tolerance must be supplied. A ridge is used only when explicitly positive.

Value

An hd_classifier_fit; non-negative scores select class 1.

References

Li, J. and Shao, J. (2015). Sparse quadratic discriminant analysis for high dimensional data. Statistica Sinica, 25, 457-473.

Examples

x <- rbind(
  c(-2, 0), c(-1, 1), c(-1, -1), c(-2, 1),
  c(2, 0), c(1, 2), c(1, -2), c(2, 1)
)
y <- factor(rep(c("left", "right"), each = 4))
fit <- li_shao_sparse_qda(
  x, y, threshold_mean = 0, threshold_difference = 0,
  threshold_covariance = 0, ridge = 0.1
)

Li–Wang–Zou simpler spatial-sign-based two-sample test

Description

Tests equality of two multivariate location vectors using the simplified bias-corrected spatial-sign test (SST) of Li, Wang, and Zou (2016). For group k=1,2, let (\widehat\theta_k,\widehat D_k) be the full-sample diagonal Hettmansperger–Randles fit satisfying

n_k^{-1}\sum_i U\{\widehat D_k^{-1/2} (X_{ki}-\widehat\theta_k)\}=0

and

p n_k^{-1}\mathrm{diag}\left(\sum_i U_{ki}U_{ki}^{\mathsf T}\right)=I_p.

Unlike the earlier Feng–Zou–Wang statistic, this method uses no leave-out fit. Its uncorrected statistic is

T_n=-\frac{1}{n_1n_2}\sum_{i=1}^{n_1}\sum_{j=1}^{n_2} U\{\widehat D_1^{-1/2}(X_{1i}-\widehat\theta_2)\}^{\mathsf T} U\{\widehat D_2^{-1/2}(X_{2j}-\widehat\theta_1)\}.

Usage

li_wang_zou_two_sample_sign_test(
  x,
  y,
  alpha = 0.05,
  tol = 1e-07,
  max_iter = 500L,
  zero_tol = 0,
  strict = TRUE
)

Arguments

Details

Write \widetilde U_{ki}=U\{\widehat D_k^{-1/2} (X_{ki}-\widehat\theta_k)\}, \widehat c_k=n_k^{-1}\sum_i \|\widehat D_k^{-1/2}(X_{ki}-\widehat\theta_k)\|^{-1}, \widehat c_{nk}=\widehat c_{3-k}/\widehat c_k, and \widehat D_n=\widehat D_1^{1/2}\widehat D_2^{1/2}. The feasible bias terms in Proposition 1 are

\widehat{\operatorname{tr}(A_k)}= \widehat c_{nk}\operatorname{tr} (\widehat D_k\widehat D_n^{-1}),\qquad \widehat\mu_n=\frac{\widehat{\operatorname{tr}(A_1)}}{n_1p}+ \frac{\widehat{\operatorname{tr}(A_2)}}{n_2p}.

The within-group squared-trace estimators are

\widehat{\operatorname{tr}(A_k^2)}= \frac{p^2\widehat c_{nk}^2}{n_k(n_k-1)} \sum_{a=1}^{n_k}\sum_{b\ne a} \{\widetilde U_{kb}^{\mathsf T}\widehat D_k \widehat D_n^{-1}\widetilde U_{ka}\}^2,

and the cross trace is

\widehat{\operatorname{tr}(A_3^{\mathsf T}A_3)}= \frac{p^2}{n_1n_2}\sum_{i=1}^{n_1}\sum_{j=1}^{n_2} (\widetilde U_{1i}^{\mathsf T}\widetilde U_{2j})^2.

The published feasible variance is

\widehat\sigma_n^2= \frac{2\widehat{\operatorname{tr}(A_1^2)}}{n_1^2p^2}+ \frac{2\widehat{\operatorname{tr}(A_2^2)}}{n_2^2p^2}+ \frac{4\widehat{\operatorname{tr}(A_3^{\mathsf T}A_3)}} {n_1n_2p^2}.

In particular, the first two outer denominators are n_k^2; only the ordered-pair trace estimators use n_k(n_k-1). The reported statistic is Z=(T_n-\widehat\mu_n)/\widehat\sigma_n.

The scientific alternative is two-sided location inequality, but the SST is quadratic in the location difference. Consequently, the original article rejects for large positive Z, and this function reports the upper standard-normal tail.

Each full-sample fit starts at its sample mean and marginal sample variances and iterates the two published estimating equations. The diagonal equations identify scale only up to a common multiplier; a unit geometric mean is imposed internally. All bias and variance combinations are invariant to that identification. iteration.stable means that both estimating-equation residuals are at most tol. With strict = TRUE, an unstable fit is an error; otherwise the last iterate is returned with a warning and full diagnostics.

The spatial-sign convention is U(0)=0 for numerator terms. A zero training radius is nevertheless an error because the recursive location update and \widehat c_k require inverse radii. No radius perturbation, weight cap, ridge, pseudoinverse, absolute-value repair, numerical floor, or variance substitution is applied.

Value

An object of class c("hd_location_test", "htest"). Its components field contains the uncorrected and bias-corrected statistics, bias and variance estimates, radial moments, three trace estimates, ordered/cross-pair kernels, both full-sample fits, signs, radii, and scale bridges. diagnostics records convergence, denominator conventions, upper-tail calibration, zero-sign use, numerical scaling, and the no-repair contract.

References

Li, Y., Wang, Z., and Zou, C. (2016). A simpler spatial-sign-based two-sample test for high-dimensional data. Journal of Multivariate Analysis, 149, 192–198. doi:10.1016/j.jmva.2016.04.004.

Examples

set.seed(2016)
x <- matrix(stats::rt(32, df = 5), 8, 4)
y <- matrix(stats::rt(36, df = 5), 9, 4)
li_wang_zou_two_sample_sign_test(x, y, tol = 1e-6)

Linear or convex pooling of covariance matrices

Description

Estimates every group covariance as a non-negative linear combination of all unbiased group SCMs and, optionally, the identity. For target group k, coefficients solve

\min_a\;\frac12a^\mathsf{T}(C+\Delta)a-C_{\cdot k}^\mathsf{T}a,

subject to the requested lower bounds. method = "convex" additionally imposes 1^\mathsf{T}a=1.

Usage

linear_pool_covariance(
  data,
  method = c("linear", "convex"),
  identity = TRUE,
  identity_lower = 0,
  tol = 1e-08,
  max_iter = 1000L,
  solver_tol = 1e-10,
  solver_max_iter = 10000L,
  strict = TRUE,
  keep_signs = FALSE
)

Arguments

Details

Each SCM uses divisor n_k-1. For group k, the spatial-median SSCM and positive residual radii give

\widehat\gamma_{0k}=\frac{pn_k}{n_k-1} \{\operatorname{tr}(S_{sgn,k}^2)-1/n_k\},\qquad \widehat\gamma_k=\Pi_{[1,p]} (\widehat\gamma_{0k}-p\widehat\delta_k),

with the Zou et al. inverse-radial bias correction. Corrected marginal kurtosis follows ollila_raninen_shrinkage_covariance(). The real-valued diagonal MSE term is

\Delta_{kk}=p^{-1}\{\tau_{1k}\operatorname{tr}(S_k)^2+ (\tau_{1k}+\tau_{2k})p\eta_k^2\gamma_k\},

where \tau_{1k}=1/(n_k-1)+\kappa_k/n_k and \tau_{2k}=\kappa_k/n_k. The diagonal of C is \eta_k^2\gamma_k; off-diagonal entries use \operatorname{tr}(S_{sgn,k}S_{sgn,l}) \operatorname{tr}(S_k)\operatorname{tr}(S_l)/p.

The independent projected-FISTA solver adds no ridge or pseudoinverse. It returns feasibility, projected-gradient, KKT, complementarity-gap, eigenvalue and restart diagnostics. Singular convex objectives are allowed. Nonconvergence is an error by default; strict = FALSE returns the last feasible iterate with a warning and labels it uncalibrated.

Value

A linear_pool_covariance list with coefficient matrix, pooled estimates, SCMs, SSCMs, centers, radii, eta/gamma/kappa, C, Delta, scaled objective matrices, and per-target solver diagnostics.

References

Raninen, E., Tyler, D. E., and Ollila, E. (2022). Linear pooling of sample covariance matrices. IEEE Transactions on Signal Processing, 70, 659–672. doi:10.1109/TSP.2021.3139207.

Examples

x1 <- matrix(c(-2, 0, -1, 2, 0, -2, 1, 1, 2, -1, 3, 2,
               -1, -2, 2, 0), ncol = 2, byrow = TRUE)
x2 <- matrix(c(-1, 1, 0, 3, 1, -2, 2, 2, 3, 0, -2, -1,
               1, 2, 0, -3), ncol = 2, byrow = TRUE)
linear_pool_covariance(list(first = x1, second = x2))

Liu–Feng–Ma robust spatial-sign sum alpha test

Description

Implements the primary heavy-tailed alpha test. Restricted factor residuals are standardized by the full-sample diagonal spatial-sign fixed point. If U_t denotes the resulting sign and h is the residualized intercept, the quadratic form is

Q=N(h'h)^{-1}\sum_{s\ne t}h_s h_t U_s'U_t.

The variance component is the primary split leave-two-out estimator. After deleting a pair, the remaining rows are split into chronological first and second halves, separate restricted slopes are fitted, and the two deleted residuals use the corresponding slopes. Its denominator is (h'h)(h'h-1), not (h'h)^2.

Usage

liu_feng_ma_spatial_sign_alpha_test(
  returns,
  factors = NULL,
  bias = c("wild_bootstrap", "none", "supplied"),
  delta_q = NULL,
  bootstrap_reps = 100L,
  seed = NULL,
  keep_bootstrap = FALSE,
  tol = 1e-07,
  max_iter = 500L,
  zero_tol = 0
)

Arguments

Details

The recommended feasible statistic is

(Q-\widehat\delta_Q)/ \{2\widehat{\operatorname{tr}(R^2)}\}^{1/2}.

By default, \widehat\delta_Q is the mean of bootstrap Q values formed by multiplying every unrestricted OLS residual entry by an independent asset-time Rademacher sign and adding the restricted fitted values. This is the method's inferential bias calibration, not a reproduction of its simulation study. An explicit seed is locally isolated: the caller's random-number state is restored even after an error. With a NULL seed the current stream advances.

The book draft omits the primary wild-bootstrap bias term, uses (h'h)^2 in the trace denominator, and describes pair-deleted scale fits rather than the paper's full-sample diagonal scale plus split pair-deleted slope fits. These choices are not silently mixed here. The book's weighted/INST construction splices this incompatible trace into a paper whose complete feasible calibration is not given in the cited material, so weighted/INST is intentionally review-only. Likewise, the dependent and Lq paragraphs do not specify a complete feasible centering, long-run variance, and calibration and are not reconstructed.

Value

An asymptotic upper-tail normal alpha-test object with every quadratic, trace, bias, scale, convergence, and bootstrap diagnostic.

References

Liu, B., Feng, L., and Ma, Y. (2023). High-dimensional alpha test of linear factor pricing models with heavy-tailed distributions. Statistica Sinica, 33, 1389–1410. doi:10.5705/ss.202021.0134.

Examples

f <- cbind(seq(-1, 1, length.out = 10))
y <- outer(seq_len(10), 1:3, function(i, j) sin(i + j / 2))
liu_feng_ma_spatial_sign_alpha_test(
  y, f, bootstrap_reps = 2, seed = 1
)

Lloyd's deterministic K-means algorithm

Description

Minimize the Chapter 7 within-cluster sum of squared Euclidean distances by alternating exact sample-mean and nearest-center updates. Exact assignment ties use the smallest cluster index. The default empty_action = "error" applies no repair. The optional "farthest" rule moves the smallest-row maximizer of current within-cluster squared distance from a donor of size at least two, then recomputes centers.

Usage

lloyd_kmeans(
  x,
  clusters,
  initial = NULL,
  initialization = c("maxmin", "random"),
  first_index = 1L,
  seed = NULL,
  ties = "first",
  empty_action = c("error", "farthest"),
  solver_tol = 1e-10,
  solver_max_iter = 100L,
  strict = TRUE,
  keep_distances = FALSE
)

Arguments

Details

initial may be a clusters by p center matrix or clusters distinct row indices. Otherwise initialization is either deterministic max-min from the explicit first_index, or seeded sampling without replacement. Seeded initialization restores the caller's RNG state.

Value

A lloyd_kmeans_fit object with labels, centers, sizes, objective trace, initialization record, and convergence/repair certificates.

References

Feng, L. (2026). High-Dimensional Data Analysis for Elliptical Symmetric Distributions, Chapter 7 (book manuscript). This implementation follows the manuscript's explicit deterministic Lloyd-update contract.

For a foundational k-means formulation, see MacQueen, J. B. (1967). Some methods for classification and analysis of multivariate observations. Proceedings of the Fifth Berkeley Symposium, 1, 281–297.

Examples

x <- rbind(c(-2, 0), c(-1, 0), c(2, 0), c(3, 0))
lloyd_kmeans(x, clusters = 2, initial = c(1, 3))

Linear programming discriminant classifier

Description

Fits the Cai–Liu linear programming discriminant (LPD) direction

\widehat\gamma\in\arg\min_\gamma\|\gamma\|_1 :\|A\gamma-(\bar x_1-\bar x_2)\|_\infty\leq\lambda,

where A is the unbiased pooled within-class covariance plus the explicitly supplied ridge I. The package score is

(z-(\bar x_1+\bar x_2)/2)^\top\widehat\gamma+ \log(\pi_1/\pi_2),

so a positive score selects class 1. The Dantzig solution is returned only after primal feasibility, dual feasibility, l1 stationarity, and a primal–dual gap are all certified.

Usage

lpd_classifier(
  x,
  y,
  lambda,
  prior = NULL,
  ridge = 0,
  solver_tol = 1e-07,
  solver_max_iter = 100000L,
  strict = TRUE,
  tie = c("class1", "class2")
)

Arguments

Value

An object of class hd_classifier_fit.

References

Cai, T. and Liu, W. (2011). A direct estimation approach to sparse linear discriminant analysis. Journal of the American Statistical Association, 106, 1566–1577. doi:10.1198/jasa.2011.tm11199.

Examples

x <- rbind(
  c(2, 1), c(1, 2), c(2, 2), c(3, 1),
  c(-2, -1), c(-1, -2), c(-2, -2), c(-3, -1)
)
y <- factor(rep(c("first", "second"), each = 4))
fit <- lpd_classifier(x, y, lambda = 0.2, ridge = 0.1)
predict(fit, x, type = "score")

Ma–Feng–Wang–Bao conditional maximum and adaptive alpha test

Description

For a supplied restricted sieve fit, the primary marginal statistics are

t_i^2=T^{-1}\widehat\sigma_{ii}^{-1} (\widehat e_{i\cdot}'1_T)^2,\qquad \widehat\sigma_{ij}=\widehat e_{i\cdot}' \widehat e_{j\cdot}/(T-d-1),

where d is the observed factor count, not the number of sieve columns. The maximum is centered by 2\log N-\log\log N and calibrated with cdf F(x)=\exp\{-\pi^{-1/2}\exp(-x/2)\}.

Usage

ma_feng_wang_bao_conditional_alpha_test(
  fit,
  factor_count = NULL,
  component = c("max", "sum", "adaptive")
)

Arguments

Details

The adaptive test combines this maximum p-value and the complete ma_lan_su_tsai_conditional_alpha_sum_test() p-value by the primary Fisher statistic -2(\log p_M+\log p_S) with a chi-squared distribution on four degrees of freedom. The book draft, and a later review paragraph, incorrectly attribute a Cauchy combination to this paper.

Value

A conditional_alpha_test/htest object retaining both component statistics, log p-values, the primary marginal divisor, and calibration diagnostics.

References

Ma, H., Feng, L., Wang, Z. and Bao, J. (2024). Adaptive testing for alphas in conditional factor models with high dimensional assets. Journal of Business & Economic Statistics, 42, 1356–1366. doi:10.1080/07350015.2024.2313543. Preprint: https://arxiv.org/abs/2307.09397.

Examples

tt <- seq(0, 1, length.out = 18)
f <- cbind(market = sin(1:18 / 3))
z <- conditional_alpha_sieve_design(cbind(tt, tt^2), f)
y <- cbind(sin(1:18), cos(1:18), sin(1:18 / 2))
fit <- conditional_alpha_sieve_fit(y, z)
ma_feng_wang_bao_conditional_alpha_test(fit, component = "max")

Ma–Lan–Su–Tsai conditional high-dimensional alpha sum test

Description

Implements the complete feasible HDA statistic for a reusable restricted sieve fit. With h=M_Z1_T and restricted residuals \widehat e_{it},

S_{NT}=\frac1{NT}\sum_i (\widehat e_{i\cdot}'1_T)^2,\qquad \widehat\mu_{NT}=\frac1{NT}\sum_{i,t} \widehat e_{it}^2h_t^2.

If q=\operatorname{ncol}(Z) and \widehat\Sigma=T^{-1}\sum_t (\widehat e_t-\bar e)(\widehat e_t-\bar e)', the implemented primary trace correction is

\widehat{\operatorname{tr}(\Sigma^2)}= \frac{T^2}{(T+q-1)(T-q)} \left\{\operatorname{tr}(\widehat\Sigma^2)- \frac{\operatorname{tr}^2(\widehat\Sigma)}{T-q}\right\}.

The variance estimate is

\widehat\sigma_{NT}^2= \frac{2\widehat{\operatorname{tr}(\Sigma^2)}}{N^2T^2} \sum_{t\ne s}h_t^2h_s^2.

Usage

ma_lan_su_tsai_conditional_alpha_sum_test(fit)

Arguments

Details

The book draft's expression using only oracle N^{-1}\operatorname{tr}(\Omega) and 2N^{-2}\operatorname{tr}(\Omega^2) omits h and the feasible finite-sample trace correction. This function uses the primary estimator and rejects non-positive estimates; it never takes an absolute value or applies a floor.

Value

A conditional_alpha_test/htest object with the upper-tail normal calibration and every feasible centering/variance component.

References

Ma, S., Lan, W., Su, L. and Tsai, C.-L. (2020). Testing alphas in conditional time-varying factor models with high-dimensional assets. Journal of Business & Economic Statistics, 38, 214–227. doi:10.1080/07350015.2018.1482758.

Examples

tt <- seq(0, 1, length.out = 18)
z <- conditional_alpha_sieve_design(cbind(tt, tt^2))
y <- cbind(sin(1:18), cos(1:18), sin(1:18 / 2))
fit <- conditional_alpha_sieve_fit(y, z)
ma_lan_su_tsai_conditional_alpha_sum_test(fit)

Mauchly–Box likelihood-ratio test of Gaussian sphericity

Description

With residual Wishart degrees of freedom m, this function computes

V=|S|/(\operatorname{tr}(S)/p)^p

and the Box–Bartlett statistic

-\left\{m-\frac{2p^2+p+2}{6p}\right\}\log V.

The denominator in the correction is p, not p+1. The latter appears in the accompanying book draft and is a transcription error.

Usage

mauchly_sphericity_test(x, center = TRUE)

Arguments

Value

An hd_covariance_test object.

References

Mauchly, J. W. (1940). Annals of Mathematical Statistics, 11, 204–209. Wang, Q. and Yao, J. (2013). Electronic Journal of Statistics, 7, 2164–2192.

Examples

set.seed(35)
x <- matrix(rnorm(60), nrow = 20, ncol = 3)
mauchly_sphericity_test(x)

Nagao's classical Gaussian identity-covariance test

Description

Tests H_0:\Sigma=I_p with

\frac{m}{2}\operatorname{tr}(S-I_p)^2,

using the covariance with divisor equal to the residual Wishart degrees of freedom m. Its fixed-dimensional reference has p(p+1)/2 degrees of freedom. The factor one-half is missing from the book draft.

Usage

nagao_identity_test(x, center = TRUE)

Arguments

Value

An hd_covariance_test object.

References

Nagao, H. (1973). On some test criteria for covariance matrix. Annals of Statistics, 1, 700–709.

Examples

set.seed(37)
x <- matrix(rnorm(60), nrow = 20, ncol = 3)
nagao_identity_test(x)

Normalize a shape matrix

Description

Removes the unidentified scalar from a positive definite shape matrix.

Usage

normalize_shape(shape, method = c("trace", "determinant"))

Arguments

Value

The symmetrized, normalized matrix.

References

Fang, K.-T. and Anderson, T. W. (1990). Statistical Inference in Elliptically Contoured and Related Distributions. Allerton Press.

Examples

shape <- matrix(c(4, 1, 1, 1), 2, 2)
normalize_shape(shape, method = "determinant")

Ollila–Raninen elliptical shrinkage covariance estimator

Description

Estimates a covariance matrix by

\widehat\Sigma=\widehat\beta S+ (1-\widehat\beta)\widehat\eta I_p,

where S is the unknown-mean unbiased sample covariance (divisor n-1), \widehat\eta=\operatorname{tr}(S)/p, and

\widehat\beta=\frac{\widehat\gamma-1} {\widehat\gamma-1+\widehat\kappa(2\widehat\gamma+p)/n+ (\widehat\gamma+p)/(n-1)}

.

Usage

ollila_raninen_shrinkage_covariance(
  x,
  sphericity = c("ell1", "ell2", "ell3"),
  tol = 1e-08,
  max_iter = 1000L,
  strict = TRUE,
  keep_signs = FALSE
)

Arguments

Details

The corrected marginal excess kurtoses are

K_j=\frac{n-1}{(n-2)(n-3)}\{(n+1)g_{2j}+6\},

and \widehat\kappa=\max\{-2/(p+2),p^{-1} \sum_jK_j/3\}. Equality at the theoretical lower bound is retained and diagnosed; the 0.99 boundary modification in the authors' MATLAB code is deliberately not used.

sphericity = "ell1" uses the spatial-median SSCM estimator

\widehat\gamma_1^*=\frac{n}{n-1} \{p\operatorname{tr}(S_{sgn}^2)-p/n\}.

"ell2" uses

\widehat\gamma_2^*=b_n\left\{ \frac{p\operatorname{tr}(S^2)}{\operatorname{tr}(S)^2} -a_n\frac pn\right\},

where a_n=n(n/(n-1)+\widehat\kappa)/(n+\widehat\kappa) and b_n=(n+\widehat\kappa)(n-1)^2/ [(n-2)\{3\widehat\kappa(n-1)+n(n+1)\}]. Both are projected to ⁠[1, p]⁠ as prescribed in the publication. "ell3" selects the smaller projected sphericity, with ties assigned to Ell2.

Value

An hd_covariance_estimator containing the estimate, unbiased SCM, scale, both raw/projected sphericities, corrected kurtosis, data weight, shrinkage intensity, and no-repair diagnostics.

References

Ollila, E. and Raninen, E. (2019). Optimal shrinkage covariance matrix estimation under random sampling from elliptical distributions. IEEE Transactions on Signal Processing, 67, 2707–2719. doi:10.1109/TSP.2019.2908144.

Examples

x <- matrix(c(-2, 0, 1, -1, 2, 0, 0, -2, 1, 1, 1, -1,
              2, 1, 0, 3, -1, 2, -2, -1, -1, 1, 3, 1),
            ncol = 3, byrow = TRUE)
ollila_raninen_shrinkage_covariance(x, sphericity = "ell3")

Oracle generic weighted-sign sum statistic

Description

Evaluates the Chapter 2 oracle statistic from a supplied reference location and diagonal shape. With

r_i=\|D^{-1/2}(X_i-\theta)\|,\quad V_i(K)=K(r_i)U\{D^{-1/2}(X_i-\theta)\},

the raw score is

T_n(K)=\frac{2}{n(n-1)}\sum_{i<j}V_i(K)^TV_j(K).

Usage

oracle_weighted_sign_sum_test(
  x,
  theta,
  diagonal,
  K = "constant",
  power = 0,
  null_sd = NULL,
  nu2 = NULL,
  trace_R2 = NULL,
  zero_tol = 0,
  keep_scores = FALSE
)

Arguments

Details

This function deliberately says oracle. It does not estimate the reference location, diagonal, radial moment, or correlation trace. A valid upper-tail normal p-value is returned only when the caller supplies either null_sd, or both nu2 and trace_R2, in which case

\sigma_{n,K}=\left\{ \frac{2\nu_{2,K}^2\operatorname{tr}(R^2)} {n(n-1)p^2}\right\}^{1/2}.

With no calibration inputs the raw score remains available and calibrated is FALSE; no feasible leave-out or plug-in estimator is guessed.

Value

A list of class oracle_weighted_sign_sum_test. The p.value field is absent unless explicit calibration is complete.

References

Feng, L., Liu, B. and Ma, Y. (2021). An inverse norm sign test for location parameters in high-dimensional data. Journal of Business & Economic Statistics 39, 807–815.

Examples

x <- matrix(c(-2, 1, 0, 3, -1, 2, 1, -3, 2, 0, 4, -2), ncol = 2)
oracle_weighted_sign_sum_test(
  x, theta = c(0, 0), diagonal = c(1, 1), K = "constant"
)

Park–Ayyala one-sample high-dimensional mean test

Description

Tests H_0: \mu=\mu_0 with the diagonal leave-two-out statistic of Park and Ayyala (2013). For every ordered pair of observations, the diagonal covariance used in the quadratic product is estimated after removing that pair. The estimated null variance uses the matching leave-two-out means. At least six observations are required: below this threshold the inverse-chi-square bias correction is not valid (and at n=5 the corrected numerator is identically zero). The algebraic API therefore permits n=6,7, but squared inverse leave-out variances have no finite Gaussian expectation at those sizes; their asymptotic normal calibration is especially fragile, so n\geq 8 is strongly preferred.

Usage

park_ayyala_one_sample_test(x, mu = NULL)

Arguments

Details

Each residual x - mu is formed in extended precision and divided by a common scale within its variable. This is algebraically neutral for the Park–Ayyala statistic and protects its leave-out moments from avoidable overflow, underflow, and cancellation under a large common translation. No zero variance is replaced, no ridge is added, and a non-positive estimated null variance is reported as an error.

Value

An object of class c("hd_location_test", "htest"). The components field contains the two ordered-pair sums, bias-corrected raw statistic, and estimated null variance. finite.sample.correction is (n-5)/(n-3), whereas bias.factor is the full coefficient (n-5)/\{n(n-1)(n-3)\} multiplying the ordered-pair sum.

References

Park, J. and Ayyala, D. N. (2013). A test for the mean vector in large dimension and small samples. Journal of Statistical Planning and Inference, 143, 929–943.

Examples

set.seed(23)
x <- matrix(rnorm(24), nrow = 8, ncol = 3)
park_ayyala_one_sample_test(x)

Pesaran CD test for cross-sectional independence

Description

Given a T by N residual matrix, computes

CD=\sqrt{2T/\{N(N-1)\}}\sum_{i<j}\hat\rho_{ij}.

The standard-normal calibration is two-sided. Residuals are used exactly as supplied: include an intercept in the preceding unit regressions when centering is required by the model.

Usage

pesaran_cd_test(residuals, keep_correlations = FALSE)

Arguments

Value

An object inheriting from htest.

References

Pesaran, M. H. (2004). General Diagnostic Tests for Cross Section Dependence in Panels. IZA Discussion Paper 1240. https://docs.iza.org/dp1240.pdf

Examples

e <- matrix(c(-2, 1, 0, 2, -1, 3, 1, -2, 2, 1, -3, 1), 4, 3)
pesaran_cd_test(e)

Pesaran–Yamagata large-N sum alpha test

Description

Implements the feasible statistic of Pesaran and Yamagata. Let v=T-K-1 and t_i be the usual unrestricted OLS intercept t statistic. Residual correlations are retained only when |\sqrt v\,\hat\rho_{ij}|> \Phi^{-1}\{1-p_0/(2N^\delta)\}. If \widetilde\rho^2 is twice the sum of squared retained upper-triangular correlations divided by N(N-1), the primary statistic is

\frac{N^{-1/2}\sum_i\{t_i^2-v/(v-2)\}} {[v/(v-2)]\{2(v-1)(1+(N-1)\widetilde\rho^2)/(v-4)\}^{1/2}}.

Usage

pesaran_yamagata_alpha_test(returns, factors = NULL, p0 = 0.1, delta = 1)

Arguments

Details

The generic dense statistic in the book draft is not this feasible finite-v test: it omits the t-square centering and the multiple-testing threshold estimator. This function follows the primary paper.

Value

An asymptotic upper-tail normal alpha-test object.

References

Pesaran, M. H. and Yamagata, T. (2023). Testing for alpha in linear factor pricing models with a large number of securities. Journal of Financial Econometrics. doi:10.1093/jjfinec/nbad002.

Examples

f <- cbind(seq(-1, 1, length.out = 12))
y <- outer(seq_len(12), 1:3, function(i, j) sin(i + j)) +
  f %*% matrix(1:3, nrow = 1)
pesaran_yamagata_alpha_test(y, f)

Penalized-matrix-decomposition sparse CCA

Description

Solves the PMD sparse CCA problem

\max_{u,v}u'Mv,\quad \lVert u\rVert_2,\lVert v\rVert_2\leq1, \quad \lVert u\rVert_1\leq c_x,\quad \lVert v\rVert_1\leq c_y

with the cross-covariance operator applied matrix-free. Later pairs use ⁠M <- M - d u v'⁠. Both actual l1 bounds must be supplied; no permutation tuning or default grid is generated.

Usage

pmd_sparse_cca(
  x,
  y,
  l1_x,
  l1_y,
  components = 1L,
  center = TRUE,
  scale = TRUE,
  covariance_divisor = c("n", "n-1"),
  initial_x = NULL,
  initial_y = NULL,
  solver_tol = 1e-08,
  solver_max_iter = 1000L,
  strict = TRUE,
  keep_operator = TRUE
)

Arguments

Details

For the shared hd_pca_fit contract, loadings, scores, center, p, and variable.names refer to the x view. The y-view counterparts are returned separately.

Value

A pmd_sparse_cca_fit inheriting from hd_pca_fit.

References

Witten, D. M., Tibshirani, R., and Hastie, T. (2009). A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis. Biostatistics, 10, 515–534.

Examples

x <- cbind(a = c(-2, -1, 1, 2), b = c(1, -1, -1, 1))
y <- cbind(c = x[, 1], d = c(-1, 1, 1, -1))
pmd_sparse_cca(x, y, l1_x = 1, l1_y = 1, scale = FALSE)

Penalized-matrix-decomposition sparse PCA

Description

For each residual data matrix R, solves the PMD rank-one problem

\max_{u,v} u'Rv,\quad \lVert u\rVert_2\leq1,\quad \lVert v\rVert_2\leq1,\quad \lVert v\rVert_1\leq c,

then applies ⁠R <- R - d u v'⁠. The supplied l1_bound is the actual c\in[1,\sqrt p], not a rescaled UI fraction. The alternating solver reports both block-KKT and fixed-point residuals.

Usage

pmd_sparse_pca(
  x,
  l1_bound,
  components = 1L,
  center = c("mean", "none"),
  scale = FALSE,
  covariance_divisor = c("n", "n-1"),
  initial = NULL,
  solver_tol = 1e-08,
  solver_max_iter = 1000L,
  symmetry_tol = sqrt(.Machine$double.eps),
  strict = TRUE,
  keep_operator = TRUE
)

Arguments

Value

A pmd_sparse_pca_fit inheriting from hd_pca_fit.

References

Witten, D. M., Tibshirani, R., and Hastie, T. (2009). A penalized matrix decomposition, with applications to sparse principal components and canonical correlation analysis. Biostatistics, 10, 515–534.

Examples

x <- rbind(c(3, 0), c(-3, 0), c(0, 1), c(0, -1))
pmd_sparse_pca(x, l1_bound = 1)

POET covariance estimator with a supplied factor count

Description

Implements Principal Orthogonal complEment Thresholding (POET). The first factors sample principal components form the low-rank part; a generalized threshold is applied only to the off-diagonal entries of the orthogonal residual covariance, and the two parts are added back. The factor count is deliberately required: the cited display in the book does not specify a unique data-driven selector.

Usage

poet_covariance(
  x,
  factors,
  threshold,
  threshold_type = c("correlation", "variance"),
  rule = c("hard", "soft", "scad", "adaptive_lasso"),
  scad_a = 3.7,
  adaptive_eta = 1,
  center = TRUE
)

Arguments

Details

threshold_type = "correlation" uses primary equation (2.6), \lambda_{ij}=\tau\sqrt{r_{ii}r_{jj}}. "variance" uses equation (3.2), C\sqrt{\hat\theta_{ij}} \{p^{-1/2}+\sqrt{\log(p)/n}\}. Residual diagonal entries are always retained. No positive-definiteness repair, ridge, or automatic tuning is applied.

Value

An hd_covariance_estimator object containing the low-rank part, raw and thresholded orthogonal complements, eigenvalues/eigenvectors, and thresholds.

References

Fan, J., Liao, Y. and Mincheva, M. (2013). Journal of the Royal Statistical Society: Series B, 75, 603–680. arXiv:1201.0175.

Examples

set.seed(34)
x <- matrix(rnorm(60), nrow = 15, ncol = 4)
poet_covariance(x, factors = 1, threshold = 0.2)

One-step Tyler POET estimator for an elliptical factor model

Description

Implements the primary POET-TME construction. With a positive-definite preliminary inverse shape \widehat V_S, it first computes

\widehat\Sigma_T=\frac{p}{n}\sum_i \frac{(X_i-\widehat\mu)(X_i-\widehat\mu)^T} {(X_i-\widehat\mu)^T\widehat V_S(X_i-\widehat\mu)}

and then applies the same supplied-factor POET decomposition once. This is the one-step estimator recommended in the primary paper, not a new convergence iteration.

Usage

poet_tme(
  x,
  factors,
  threshold = NULL,
  constant = 1,
  rule = c("hard", "soft", "scad", "adaptive_lasso"),
  scad_a = 3.7,
  adaptive_eta = 1,
  preliminary_precision = NULL,
  tol = 1e-08,
  max_iter = 500L,
  zero_tol = 0,
  strict = TRUE
)

Arguments