Pn: Utilities for projected-ecdf statistics of spherical uniformity

Description

Computation of the kernels $$\psi_p^W(\theta) := \int_{-1}^1 A_x(\theta)\,\mathrm{d}W(F_p(x)),$$ where $A_x(\theta)$ is the proportion of area surface of $S^{p - 1}$ covered by the intersection of two hyperspherical caps with common solid angle $\pi - \cos^{-1}(x)$ and centers separated by an angle $\theta \in [0, \pi]$, $F_p$ is the distribution function of the projected spherical uniform distribution, and $W$ is a measure on $[0, 1]$.

Also, computation of the Gegenbauer coefficients of $\psi_p^W$: $$b_{k, p}^W := \frac{1}{c_{k, p}}\int_0^\pi \psi_p^W(\theta) C_k^{p / 2 - 1}(\cos\theta)\,\mathrm{d}\theta.$$ These coefficients can also be computed via $$b_{k, p}^W = \int_{-1}^1 a_{k, p}^x\,\mathrm{d}W(F_p(x))$$ for a certain function $x\rightarrow a_{k, p}^x$. They serve to define projected alternatives to uniformity.

Usage

psi_Pn(theta, q, type, Rothman_t = 1/3, tilde = FALSE, psi_Gauss = TRUE,
  psi_N = 320, tol = 1e-06)
Gegen_coefs_Pn(k, p, type, Rothman_t = 1/3, Gauss = TRUE, N = 320,
  tol = 1e-06, verbose = FALSE)
akx(x, p, k, sqr = FALSE)
f_locdev_Pn(p, type, K = 1000, N = 320, K_max = 10000, thre = 0.001,
  Rothman_t = 1/3, verbose = FALSE)

Value

psi_Pn: a vector of size length(theta) with the evaluation of $\psi$.
Gegen_coefs_Pn: a vector of size length(k) containing the coefficients $b_{k, p}^W$.
akx: a matrix of size c(length(x), length(k)) containing the coefficients $a_{k, p}^x$.
f_locdev_Pn: the projected alternative $f$ as a function ready to be evaluated.

Arguments

theta: vector with values in $[0, \pi]$.
q: integer giving the dimension of the sphere $S^q$.
type: type of projected-ecdf test statistic. Must be either "PCvM" (Cramér--von Mises), "PAD" (Anderson--Darling), or "PRt" (Rothman).
Rothman_t: $t$ parameter for the Rothman test, a real in $(0, 1)$. Defaults to 1 / 3.
tilde: include the constant and bias term? Defaults to FALSE.
psi_Gauss: use a Gauss--Legendre quadrature rule with psi_N nodes in the computation of the kernel function? Defaults to TRUE.
psi_N: number of points used in the Gauss--Legendre quadrature for computing the kernel function. Defaults to 320.
tol: tolerance passed to integrate's rel.tol and abs.tol if Gauss = FALSE. Defaults to 1e-6.
k: vector with the index of coefficients.
p: integer giving the dimension of the ambient space $R^p$ that contains $S^{p-1}$.
Gauss: use a Gauss--Legendre quadrature rule of N nodes in the computation of the Gegenbauer coefficients? Otherwise, call integrate. Defaults to TRUE.
N: number of points used in the Gauss--Legendre quadrature for computing the Gegenbauer coefficients. Defaults to 320.
verbose: flag to print informative messages. Defaults to FALSE.
x: evaluation points for $a_{k, p}^x$, a vector with values in $[-1, 1]$.
sqr: return the signed square root of $a_{k, p}^x$? Defaults to FALSE.
K: number of equispaced points on $[-1, 1]$ used for evaluating $f$ and then interpolating. Defaults to 1e3.
K_max: integer giving the truncation of the series. Defaults to 1e4.
thre: proportion of norm not explained by the first terms of the truncated series. Defaults to 1e-3.

Author

Eduardo García-Portugués and Paula Navarro-Esteban.

Details

The evaluation of $\psi_p^W$ and $b_{k, p}^W$ depends on the type of projected-ecdf statistic:

PCvM: closed-form expressions for $\psi_p^W$ and $b_{k, p}^W$ with $p = 2, 3, 4$, numerical integration required for $p \ge 5$.
PAD: closed-form expressions for $\psi_2^W$ and $b_{k, 3}^W$, numerical integration required for $\psi_p^W$ with $p \ge 3$ and $b_{k, p}^W$ with $p = 2$ and $p \ge 4$.
PRt: closed-form expressions for $\psi_p^W$ and $b_{k, p}^W$ for any $p \ge 2$.

See García-Portugués et al. (2023) for more details.

References

García-Portugués, E., Navarro-Esteban, P., Cuesta-Albertos, J. A. (2023) On a projection-based class of uniformity tests on the hypersphere. Bernoulli, 29(1):181--204. tools:::Rd_expr_doi("10.3150/21-BEJ1454")

Examples

Run this code

# Kernels in the projected-ecdf test statistics
k <- 0:10
coefs <- list()
(coefs$PCvM <- t(sapply(2:5, function(p)
  Gegen_coefs_Pn(k = k, p = p, type = "PCvM"))))
(coefs$PAD <- t(sapply(2:5, function(p)
  Gegen_coefs_Pn(k = k, p = p, type = "PAD"))))
(coefs$PRt <- t(sapply(2:5, function(p)
  Gegen_coefs_Pn(k = k, p = p, type = "PRt"))))

# Gegenbauer expansion
th <- seq(0, pi, length.out = 501)[-501]
old_par <- par(mfrow = c(3, 4))
for (type in c("PCvM", "PAD", "PRt")) {

  for (p in 2:5) {

    plot(th, psi_Pn(theta = th, q = p - 1, type = type), type = "l",
         main = paste0(type, ", p = ", p), xlab = expression(theta),
         ylab = expression(psi(theta)), axes = FALSE, ylim = c(-1.5, 1))
    axis(1, at = c(0, pi / 4, pi / 2, 3 * pi / 4, pi),
         labels = expression(0, pi / 4, pi / 2, 3 * pi / 4, pi))
    axis(2); box()
    lines(th, Gegen_series(theta = th, coefs = coefs[[type]][p - 1, ],
                           k = k, p = p), col = 2)

  }

}
par(old_par)

# Analytical coefficients vs. numerical integration
test_coef <- function(type, p, k = 0:20) {

  plot(k, log1p(abs(Gegen_coefs_Pn(k = k, p = p, type = type))),
       ylab = "Coefficients", main = paste0(type, ", p = ", p))
  points(k, log1p(abs(Gegen_coefs(k = k, p = p, psi = psi_Pn, type = type,
                                  q = p - 1))), col = 2)
  legend("topright", legend = c("log(1 + Gegen_coefs_Pn))",
                                "log(1 + Gegen_coefs(psi_Pn))"),
         lwd = 2, col = 1:2)

}

# PCvM statistic
old_par <- par(mfrow = c(2, 2))
for (p in 2:5) test_coef(type = "PCvM", p = p)
par(old_par)

# PAD statistic
old_par <- par(mfrow = c(2, 2))
for (p in 2:5) test_coef(type = "PAD", p = p)
par(old_par)

# PRt statistic
old_par <- par(mfrow = c(2, 2))
for (p in 2:5) test_coef(type = "PRt", p = p)
par(old_par)

# akx
akx(x = seq(-1, 1, l = 5), k = 1:4, p = 2)
akx(x = 0, k = 1:4, p = 3)

# PRt alternative to uniformity
z <- seq(-1, 1, l = 1e3)
p <- c(2:5, 10, 15, 17)
col <- viridisLite::viridis(length(p))
plot(z, f_locdev_Pn(p = p[1], type = "PRt")(z), type = "s",
     col = col[1], ylim = c(0, 0.6), ylab = expression(f[Rt](z)))
for (k in 2:length(p)) {
  lines(z, f_locdev_Pn(p = p[k], type = "PRt")(z), type = "s", col = col[k])
}
legend("topleft", legend = paste("p =", p), col = col, lwd = 2)

Run the code above in your browser using DataLab