Pn: Utilities for projected-ecdf statistics of spherical...

In sphunif: Uniformity Tests on the Circle, Sphere, and Hypersphere

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)

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.

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.

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.

Author(s)

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

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. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.3150/21-BEJ1454")}.

Examples

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

sphunif documentation built on May 29, 2024, 4:19 a.m.