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

Pn {sphunif}

R Documentation

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

[Package sphunif version 1.4.0 Index]