Package 'sphunif'

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

Description

Implementation of uniformity tests on the circle and (hyper)sphere. The main function of the package is unif_test, which conveniently collects more than 35 tests for assessing uniformity on $S^{p-1}=\{{\bf x}\in R^p:||{\bf x}||=1\}$, $p\ge 2$. The test statistics are implemented in the unif_stat function, which allows computing several statistics for different samples within a single call, thus facilitating Monte Carlo experiments. Furthermore, the unif_stat_MC function allows parallelizing them in a simple way. The asymptotic null distributions of the statistics are available through the function unif_stat_distr. The core of sphunif-package is coded in C++ by relying on the Rcpp-package. The package also provides several novel datasets and gives the replicability for the data applications/ simulations in García-Portugués et al. (2021) <doi:10.1007/978-3-030-69944-4_12>, García-Portugués et al. (2023) <doi:10.3150/21-BEJ1454>, García-Portugués et al. (2024) <doi:10.48550/arXiv.2108.09874>, and Fernández-de-Marcos and García-Portugués (2024) <doi:10.48550/arXiv.405.13531>.

Author(s)

Eduardo García-Portugués and Thomas Verdebout.

References

Fernández-de-Marcos, A. and García-Portugués, E. (2024) A stereographic test of spherical uniformity. arXiv:2405.13531. doi:10.48550/arXiv.2405.13531.

García-Portugués, E. and Verdebout, T. (2018) An overview of uniformity tests on the hypersphere. arXiv:1804.00286. doi:10.48550/arXiv.1804.00286.

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.

García-Portugués, E., Navarro-Esteban, P., and Cuesta-Albertos, J. A. (2021). A Cramér–von Mises test of uniformity on the hypersphere. In Balzano, S., Porzio, G. C., Salvatore, R., Vistocco, D., and Vichi, M. (Eds.), Statistical Learning and Modeling in Data Analysis, Studies in Classification, Data Analysis and Knowledge Organization, pp. 107–116. Springer, Cham. doi:10.1007/978-3-030-69944-4_12.

García-Portugués, E., Paindaveine, D., and Verdebout, T. (2024). On a class of Sobolev tests for symmetry of directions, their detection thresholds, and asymptotic powers. arXiv:2108.09874v2. doi:10.48550/arXiv.2108.09874.

---

Surface area of the intersection of two hyperspherical caps

Description

Computation of

$A_x(\theta_{ij}) := \frac{1}{\omega_p} \int_{S^{p - 1}} 1_{\{{\bf X}_i'\boldsymbol\gamma \le x, {\bf X}_j'\boldsymbol\gamma \le x\}}\,\mathrm{d}\boldsymbol\gamma,$

where $\theta_{ij} := \cos^{-1}({\bf X}_i'{\bf X}_j) \in [0, \pi]$, $x \in [-1, 1]$, and $\omega_{p}$ is the surface area of $S^{p - 1}$. $A_x(\theta_{ij})$ is the proportion of surface area of $S^{p - 1}$ covered by the intersection of two hyperspherical caps centered at ${\bf X}_i$ and ${\bf X}_j$ and with common solid angle $\pi - \cos^{-1}(x)$.

Usage

A_theta_x(theta, x, p, N = 160L, as_matrix = TRUE)

Arguments

theta	vector with values in $[0, \pi]$.
x	vector with values in $[-1, 1]$.
p	integer giving the dimension of the ambient space $R^p$ that contains $S^{p-1}$.
N	number of points used in the Gauss-Legendre quadrature. Defaults to 160.
as_matrix	return a matrix with the values of $A_x(\theta)$ on the grid formed by theta and x? If FALSE, $A_x(\theta)$ is evaluated on theta and x if they equal in size. Defaults to TRUE.

Details

See García-Portugués et al. (2023) for more details about the $A_x(\theta)$ function.

Value

A matrix of size c(length(theta), length(x)) containing the evaluation of $A_x(\theta)$ if as_matrix = TRUE. Otherwise, a vector of size c(length(theta) if theta and x equal in size.

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

# Plot A_x(theta) for several dimensions and x's
A_lines <- function(x, th = seq(0, pi, l = 200)) {

  plot(th, A_theta_x(theta = th, x = x, p = 2), type = "l",
       col = 1, ylim = c(0, 1.25), main = paste("x =", x),
       ylab = expression(A[x](theta)),
       xlab = expression(theta), axes = FALSE)
  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()
  abline(h = c(0, 1), lty = 2)
  lines(th, A_theta_x(theta = th, x = x, p = 3), col = 2)
  lines(th, A_theta_x(theta = th, x = x, p = 4), col = 3)
  lines(th, A_theta_x(theta = th, x = x, p = 5), col = 4)
  legend("top", lwd = 2, legend = paste("p =", 2:5),
         col = 1:4, cex = 0.75, horiz = TRUE)

}
old_par <- par(mfrow = c(2, 3))
A_lines(x = -0.75)
A_lines(x = -0.25)
A_lines(x = 0)
A_lines(x = 0.25)
A_lines(x = 0.5)
A_lines(x = 0.75)
par(old_par)

# As surface of (theta, x) for several dimensions
A_surf <- function(p, x = seq(-1, 1, l = 201), th = seq(0, pi, l = 201)) {

  col <- c("white", viridisLite::viridis(20))
  breaks <- c(-1, seq(1e-15, 1, l = 21))
  A <- A_theta_x(theta = th, x = x, p = p)
  image(th, x, A, main = paste("p =", p), col = col, breaks = breaks,
        xlab = expression(theta), axes = FALSE)
  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()
  contour(th, x, A, levels = breaks, add = TRUE)

}
old_par <- par(mfrow = c(2, 2))
A_surf(p = 2)
A_surf(p = 3)
A_surf(p = 4)
A_surf(p = 5)
par(old_par)

# No matrix return
th <- seq(0, pi, l = 5)
x <- seq(-1, 1, l = 5)
diag(A_theta_x(theta = th, x = x, p = 2))
A_theta_x(theta = th, x = x, p = 2, as_matrix = FALSE)

---

Conversion between angular and Cartesian coordinates of the (hyper)sphere

Description

Transforms the angles $(\theta_1,\ldots,\theta_{p-1})'$ in $[0,\pi)^{p-2}\times[-\pi,\pi)$ into the Cartesian coordinates

$(\cos(x_1),\sin(x_1)\cos(x_2),\ldots, \sin(x_1)\cdots\sin(x_{p-2})\cos(x_{p-1}), \sin(x_1)\cdots\sin(x_{p-2})\sin(x_{p-1}))'$

of $S^{p-1}$, and vice versa.

Usage

angles_to_sphere(theta)

sphere_to_angles(x)

Arguments

theta	matrix of size c(n, p - 1) with the angles.
x	matrix of size c(n, p) with the Cartesian coordinates. Assumed to be of unit norm by rows.

Value

For angles_to_sphere, the matrix x. For sphere_to_angles, the matrix theta.

Examples

# Check changes of coordinates
sphere_to_angles(angles_to_sphere(c(pi / 2, 0, pi)))
sphere_to_angles(angles_to_sphere(rbind(c(pi / 2, 0, pi), c(pi, pi / 2, 0))))
angles_to_sphere(sphere_to_angles(c(0, sqrt(0.5), sqrt(0.1), sqrt(0.4))))
angles_to_sphere(sphere_to_angles(
  rbind(c(0, sqrt(0.5), sqrt(0.1), sqrt(0.4)),
        c(0, sqrt(0.5), sqrt(0.5), 0),
        c(0, 1, 0, 0),
        c(0, 0, 0, -1),
        c(0, 0, 1, 0))))

# Circle
sphere_to_angles(angles_to_sphere(0))
sphere_to_angles(angles_to_sphere(cbind(0:3)))
angles_to_sphere(cbind(sphere_to_angles(rbind(c(0, 1), c(1, 0)))))
angles_to_sphere(cbind(sphere_to_angles(rbind(c(0, 1)))))

---

Available circular and (hyper)spherical uniformity tests

Description

Listing of the tests implemented in the sphunif package.

Usage

avail_cir_tests

avail_sph_tests

Format

An object of class character of length 33.

An object of class character of length 18.

Value

A character vector whose elements are valid inputs for the type argument in unif_test, unif_stat, unif_stat_distr, and unif_stat_MC. avail_cir_tests provides the available circular tests and avail_sph_tests the (hyper)spherical tests.

Examples

# Circular tests
avail_cir_tests

# Spherical tests
avail_sph_tests

---

Transforming between polar and Cartesian coordinates

Description

Transformation between a matrix Theta containing M circular samples of size n on $[0, 2\pi)$ and an array X containing the associated Cartesian coordinates on $S^1:=\{{\bf x}\in R^2:||{\bf x}||=1\}$.

Usage

Theta_to_X(Theta)

X_to_Theta(X)

Arguments

Theta	a matrix of size c(n, M) with M samples of size n of circular data on $[0, 2\pi)$. Must not contain NA's.
X	an array of size c(n, 2, M) containing the Cartesian coordinates of M samples of size n of directions on $S^{1}$. Must not contain NA's.

Value

• Theta_to_X: the corresponding X.

• X_to_Theta: the corresponding Theta.

Examples

# Sample
Theta <- r_unif_cir(n = 10, M = 2)
X <- r_unif_sph(n = 10, p = 2, M = 2)

# Check equality
sum(abs(X - Theta_to_X(X_to_Theta(X))))
sum(abs(Theta - X_to_Theta(Theta_to_X(Theta))))

---

Circular gaps

Description

Computation of the circular gaps of an angular sample $\Theta_1,\ldots,\Theta_n$ on $[0, 2\pi)$, defined as

$\Theta_{(2)} - \Theta_{(1)},\ldots,\Theta_{(n)} - \Theta_{(n - 1)}, 2\pi - \Theta_{(n)} - \Theta_{(1)},$

where

$0 \le \Theta_{(1)} \le \Theta_{(2)} \le \ldots \le \Theta_{(n)} \le 2\pi.$

Usage

cir_gaps(Theta, sorted = FALSE)

Arguments

Theta	a matrix of size c(n, M) with M samples of size n of circular data on $[0, 2\pi)$. Must not contain NA's.
sorted	are the columns of Theta sorted increasingly? If TRUE, performance is improved. If FALSE (default), each column of Theta is sorted internally.

Value

A matrix of size c(n, M) containing the n circular gaps for each of the M circular samples.

Warning

Be careful on avoiding the next bad usages of cir_gaps, which will produce spurious results:

• The entries of Theta are not in $[0, 2\pi)$.

• Theta is not sorted increasingly when data_sorted = TRUE.

Examples

Theta <- cbind(c(pi, 0, 3 * pi / 2), c(0, 3 * pi / 2, pi), c(5, 3, 1))
cir_gaps(Theta)

---

Statistics for testing circular uniformity

Description

Low-level implementation of several statistics for assessing circular uniformity on $[0, 2\pi)$ or, equivalently, $S^1:=\{{\bf x}\in R^2:||{\bf x}||=1\}$.

Usage

cir_stat_Kuiper(Theta, sorted = FALSE, KS = FALSE, Stephens = FALSE)

cir_stat_Watson(Theta, sorted = FALSE, CvM = FALSE, Stephens = FALSE)

cir_stat_Watson_1976(Theta, sorted = FALSE, minus = FALSE)

cir_stat_Range(Theta, sorted = FALSE, gaps_in_Theta = FALSE,
  max_gap = TRUE)

cir_stat_Rao(Theta, sorted = FALSE, gaps_in_Theta = FALSE)

cir_stat_Greenwood(Theta, sorted = FALSE, gaps_in_Theta = FALSE)

cir_stat_Log_gaps(Theta, sorted = FALSE, gaps_in_Theta = FALSE,
  abs_val = TRUE)

cir_stat_Vacancy(Theta, a = 2 * pi, sorted = FALSE,
  gaps_in_Theta = FALSE)

cir_stat_Max_uncover(Theta, a = 2 * pi, sorted = FALSE,
  gaps_in_Theta = FALSE)

cir_stat_Num_uncover(Theta, a = 2 * pi, sorted = FALSE,
  gaps_in_Theta = FALSE, minus_val = TRUE)

cir_stat_Gini(Theta, sorted = FALSE, gaps_in_Theta = FALSE)

cir_stat_Gini_squared(Theta, sorted = FALSE, gaps_in_Theta = FALSE)

cir_stat_Ajne(Theta, Psi_in_Theta = FALSE)

cir_stat_Rothman(Theta, Psi_in_Theta = FALSE, t = 1/3)

cir_stat_Hodges_Ajne(Theta, asymp_std = FALSE, sorted = FALSE,
  use_Cressie = TRUE)

cir_stat_Cressie(Theta, t = 1/3, sorted = FALSE)

cir_stat_FG01(Theta, sorted = FALSE)

cir_stat_Rayleigh(Theta, m = 1L)

cir_stat_Bingham(Theta)

cir_stat_Hermans_Rasson(Theta, Psi_in_Theta = FALSE)

cir_stat_Gine_Gn(Theta, Psi_in_Theta = FALSE)

cir_stat_Gine_Fn(Theta, Psi_in_Theta = FALSE)

cir_stat_Pycke(Theta, Psi_in_Theta = FALSE)

cir_stat_Pycke_q(Theta, Psi_in_Theta = FALSE, q = 0.5)

cir_stat_Bakshaev(Theta, Psi_in_Theta = FALSE)

cir_stat_Riesz(Theta, Psi_in_Theta = FALSE, s = 1)

cir_stat_PCvM(Theta, Psi_in_Theta = FALSE)

cir_stat_PRt(Theta, Psi_in_Theta = FALSE, t = 1/3)

cir_stat_PAD(Theta, Psi_in_Theta = FALSE, AD = FALSE, sorted = FALSE)

cir_stat_Poisson(Theta, Psi_in_Theta = FALSE, rho = 0.5)

cir_stat_Softmax(Theta, Psi_in_Theta = FALSE, kappa = 1)

cir_stat_CCF09(Theta, dirs, K_CCF09 = 25L, original = FALSE)

Arguments

Theta	a matrix of size c(n, M) with M samples of size n of circular data on $[0, 2\pi)$. Must not contain NA's.
sorted	are the columns of Theta sorted increasingly? If TRUE, performance is improved. If FALSE (default), each column of Theta is sorted internally.
KS	compute the Kolmogorov-Smirnov statistic (which is not invariant under origin shifts) instead of the Kuiper statistic? Defaults to FALSE.
Stephens	compute Stephens (1970) modification so that the null distribution of the is less dependent on the sample size? The modification does not alter the test decision.
CvM	compute the Cramér-von Mises statistic (which is not invariant under origin shifts) instead of the Watson statistic? Defaults to FALSE.
minus	compute the invariant $D_n^-$ instead of $D_n^+$? Defaults to FALSE.
gaps_in_Theta	does Theta contain the matrix of circular gaps that is obtained with cir_gaps(Theta)? If FALSE (default), the circular gaps are computed internally.
max_gap	compute the maximum gap for the range statistic? If TRUE (default), rejection happens for large values of the statistic, which is consistent with the rest of tests. Otherwise, the minimum gap is computed and rejection happens for low values.
abs_val	return the absolute value of the Darling's log gaps statistic? If TRUE (default), rejection happens for large values of the statistic, which is consistent with the rest of tests. Otherwise, the signed statistic is computed and rejection happens for large absolute values.
a	either: $a_n = a / n$ parameter used in the length of the arcs of the coverage-based tests. Must be positive. Defaults to 2 * pi. $a$ parameter for the Stereo test, a real in $[-1, 1]$. Defaults to 0.
minus_val	return the negative value of the (standardized) number of uncovered spacings? If TRUE (default), rejection happens for large values of the statistic, which is consistent with the rest of tests. Otherwise, rejection happens for low values.
Psi_in_Theta	does Theta contain the shortest angles matrix $\boldsymbol\Psi$ that is obtained with Psi_mat(array(Theta, dim = c(n, 1, M)))? If FALSE (default), $\boldsymbol\Psi$ is computed internally.
t	$t$ parameter for the Rothman and Cressie tests, a real in $(0, 1)$. Defaults to 1 / 3.
asymp_std	normalize the Hodges-Ajne statistic in terms of its asymptotic distribution? Defaults to FALSE.
use_Cressie	compute the Hodges-Ajne statistic as a particular case of the Cressie statistic? Defaults to TRUE as it is more efficient. If FALSE, the geometric construction in Ajne (1968) is employed.
m	integer $m$ for the $m$-modal Rayleigh test. Defaults to m = 1 (the standard Rayleigh test).
q	$q$ parameter for the Pycke "$q$-test", a real in $(0, 1)$. Defaults to 1 / 2.
s	$s$ parameter for the $s$-Riesz test, a real in $(0, 2)$. Defaults to 1.
AD	compute the Anderson-Darling statistic (which is not invariant under origin shifts) instead of the Projected Anderson-Darling statistic? Defaults to FALSE.
rho	$\rho$ parameter for the Poisson test, a real in $[0, 1)$. Defaults to 0.5.
kappa	$\kappa$ parameter for the Softmax test, a non-negative real. Defaults to 1.
dirs	a matrix of size c(n_proj, 2) containing n_proj random directions (in Cartesian coordinates) on $S^1$ to perform the CCF09 test.
K_CCF09	integer giving the truncation of the series present in the asymptotic distribution of the Kolmogorov-Smirnov statistic. Defaults to 25.
original	return the CCF09 statistic as originally defined? If FALSE (default), a faster and equivalent statistic is computed, and rejection happens for large values of the statistic, which is consistent with the rest of tests. Otherwise, rejection happens for low values.

Details

Descriptions and references for most of the statistics are available in García-Portugués and Verdebout (2018).

The statistics cir_stat_PCvM and cir_stat_PRt are provided for the sake of completion, but they equal the more efficiently-implemented statistics 2 * cir_stat_Watson and cir_stat_Rothman, respectively.

Value

A matrix of size c(M, 1) containing the statistics for each of the M samples.

Warning

Be careful on avoiding the next bad usages of the functions, which will produce spurious results:

• The entries of Theta are not in $[0, 2\pi)$.

• Theta does not contain the circular gaps when gaps_in_Theta = TRUE.

• Theta is not sorted increasingly when data_sorted = TRUE.

• Theta does not contain Psi_mat(array(Theta, dim = c(n, 1, M))) when
  Psi_in_Theta = TRUE.

• The directions in dirs do not have unit norm.

References

García-Portugués, E. and Verdebout, T. (2018) An overview of uniformity tests on the hypersphere. arXiv:1804.00286. doi:10.48550/arXiv.1804.00286.

Examples

## Sample uniform circular data

M <- 2
n <- 100
set.seed(987202226)
Theta <- r_unif_cir(n = n, M = M)

## Tests based on the empirical cumulative distribution function

# Kuiper
cir_stat_Kuiper(Theta)
cir_stat_Kuiper(Theta, Stephens = TRUE)

# Watson
cir_stat_Watson(Theta)
cir_stat_Watson(Theta, Stephens = TRUE)

# Watson (1976)
cir_stat_Watson_1976(Theta)

## Partition-based tests

# Ajne
Theta_array <- Theta
dim(Theta_array) <- c(nrow(Theta), 1, ncol(Theta))
Psi <- Psi_mat(Theta_array)
cir_stat_Ajne(Theta)
cir_stat_Ajne(Psi, Psi_in_Theta = TRUE)

# Rothman
cir_stat_Rothman(Theta, t = 0.5)
cir_stat_Rothman(Theta)
cir_stat_Rothman(Psi, Psi_in_Theta = TRUE)

# Hodges-Ajne
cir_stat_Hodges_Ajne(Theta)
cir_stat_Hodges_Ajne(Theta, use_Cressie = FALSE)

# Cressie
cir_stat_Cressie(Theta, t = 0.5)
cir_stat_Cressie(Theta)

# FG01
cir_stat_FG01(Theta)

## Spacings-based tests

# Range
cir_stat_Range(Theta)

# Rao
cir_stat_Rao(Theta)

# Greenwood
cir_stat_Greenwood(Theta)

# Log gaps
cir_stat_Log_gaps(Theta)

# Vacancy
cir_stat_Vacancy(Theta)

# Maximum uncovered spacing
cir_stat_Max_uncover(Theta)

# Number of uncovered spacings
cir_stat_Num_uncover(Theta)

# Gini mean difference
cir_stat_Gini(Theta)

# Gini mean squared difference
cir_stat_Gini_squared(Theta)

## Sobolev tests

# Rayleigh
cir_stat_Rayleigh(Theta)
cir_stat_Rayleigh(Theta, m = 2)

# Bingham
cir_stat_Bingham(Theta)

# Hermans-Rasson
cir_stat_Hermans_Rasson(Theta)
cir_stat_Hermans_Rasson(Psi, Psi_in_Theta = TRUE)

# Gine Fn
cir_stat_Gine_Fn(Theta)
cir_stat_Gine_Fn(Psi, Psi_in_Theta = TRUE)

# Gine Gn
cir_stat_Gine_Gn(Theta)
cir_stat_Gine_Gn(Psi, Psi_in_Theta = TRUE)

# Pycke
cir_stat_Pycke(Theta)
cir_stat_Pycke(Psi, Psi_in_Theta = TRUE)

# Pycke q
cir_stat_Pycke_q(Theta)
cir_stat_Pycke_q(Psi, Psi_in_Theta = TRUE)

# Bakshaev
cir_stat_Bakshaev(Theta)
cir_stat_Bakshaev(Psi, Psi_in_Theta = TRUE)

# Riesz
cir_stat_Riesz(Theta, s = 1)
cir_stat_Riesz(Psi, Psi_in_Theta = TRUE, s = 1)

# Projected Cramér-von Mises
cir_stat_PCvM(Theta)
cir_stat_PCvM(Psi, Psi_in_Theta = TRUE)

# Projected Rothman
cir_stat_PRt(Theta, t = 0.5)
cir_stat_PRt(Theta)
cir_stat_PRt(Psi, Psi_in_Theta = TRUE)

# Projected Anderson-Darling
cir_stat_PAD(Theta)
cir_stat_PAD(Psi, Psi_in_Theta = TRUE)

## Other tests

# CCF09
dirs <- r_unif_sph(n = 3, p = 2, M = 1)[, , 1]
cir_stat_CCF09(Theta, dirs = dirs)

## Connection of Kuiper and Watson statistics with KS and CvM, respectively

# Rotate sample for KS and CvM
alpha <- seq(0, 2 * pi, l = 1e4)
KS_alpha <- sapply(alpha, function(a) {
  cir_stat_Kuiper((Theta[, 2, drop = FALSE] + a) %% (2 * pi), KS = TRUE)
})
CvM_alpha <- sapply(alpha, function(a) {
  cir_stat_Watson((Theta[, 2, drop = FALSE] + a) %% (2 * pi), CvM = TRUE)
})
AD_alpha <- sapply(alpha, function(a) {
  cir_stat_PAD((Theta[, 2, drop = FALSE] + a) %% (2 * pi), AD = TRUE)
})

# Kuiper is the maximum rotated KS
plot(alpha, KS_alpha, type = "l")
abline(h = cir_stat_Kuiper(Theta[, 2, drop = FALSE]), col = 2)
points(alpha[which.max(KS_alpha)], max(KS_alpha), col = 2, pch = 16)

# Watson is the minimum rotated CvM
plot(alpha, CvM_alpha, type = "l")
abline(h = cir_stat_Watson(Theta[, 2, drop = FALSE]), col = 2)
points(alpha[which.min(CvM_alpha)], min(CvM_alpha), col = 2, pch = 16)

# Anderson-Darling is the average rotated AD?
plot(alpha, AD_alpha, type = "l")
abline(h = cir_stat_PAD(Theta[, 2, drop = FALSE]), col = 2)
abline(h = mean(AD_alpha), col = 3)

---

Comet orbits

Description

Comet orbits data from the JPL Small-Body Database Search Engine. The normal vector of a comet orbit represents is a vector on $S^2$.

Usage

comets

Format

A data frame with 3798 rows and 13 variables:

id

database ID.

spkid

object primary SPK-ID.

full_name

full name/designation following the IUA naming convention.

pdes

object primary designation.

frag

flag indicating if the record is a comet fragment.

diameter

diameter from equivalent sphere (in km).

i

inclination; the orbit's plane angle with respect to the ecliptic plane, in radians in $[0, \pi]$.

om

longitude of the ascending node; the counterclockwise angle from the vector pointing to the First Point of Aries and that pointing to the ascending node (the intersection between orbit and ecliptic plane), in radians in $[0, 2\pi)$. (Both vectors are heliocentric and within the ecliptic plane.)

per_y

sidereal orbital period (in years).

class

orbit classification. A factor with levels given below.

e

eccentricity of the orbit.

a

semi-major axis of the orbit (in AU).

w

argument of perihelion; the (shortest) angle between the vector pointing to the ascending node and that pointing to the perihelion (nearest orbit point to the Sun), in radians in $[0, \pi]$. (Both vectors are heliocentric and within the orbit's plane.)

first_obs, last_obs

Date of the first and last recorded observations used in the orbit fit.

ccf09

flag indicating if the comet was considered in the data application in Cuesta-Albertos et al. (2009); see details below.

Details

The normal vector to the ecliptic plane of the comet with inclination $i$ and longitude of the ascending node $\omega$ is

$(\sin(i) \sin(\omega), -\sin(i) \cos(\omega), \cos(i))'.$

A prograde comet has positive $\cos(i)$, negative $\cos(i)$ represents a retrograde comet.

class has the following levels:

• COM: comet orbit not matching any defined orbit class.

• CTc: Chiron-type comet, as defined by Levison and Duncan (T_Jupiter > 3; a > a_Jupiter).

• ETc: Encke-type comet, as defined by Levison and Duncan (T_Jupiter > 3; a < a_Jupiter).

• HTC: Halley-type comet, classical definition (20y < P < 200y).

• HYP: comets on hyperbolic orbits.

• JFc: Jupiter-family comet, as defined by Levison and Duncan (2 < T_Jupiter < 3).

• JFC: Jupiter-family comet, classical definition (P < 20y).

• PAR: comets on parabolic orbits.

Hyperbolic and parabolic comets are not periodic; only elliptical comets are periodic.

The ccf09 variable gives the observations considered in Cuesta-Albertos et al. (2009) after fetching in the database in 2007-12-14 for the comets such that !(class %in% c("HYP", "PAR")) & per_y >= 200. Due to the dynamic nature of the data, more comets were added to the database since 2007 and also some past records were updated.

The script performing the data preprocessing is available at comets.R. The data was retrieved on 2022-05-28. A previous version of this dataset based on the old NASA's JPL Database (accessed on 2020-05-07) is available at comets-old.rda and was obtained with comets-old.R.

Source

https://ssd.jpl.nasa.gov/tools/sbdb_query.html

References

Cuesta-Albertos, J. A., Cuevas, A., Fraiman, R. (2009) On projection-based tests for directional and compositional data. Statistics and Computing, 19:367–380. doi:10.1007/s11222-008-9098-3

Examples

# Load data
data("comets")

# Add normal vectors
comets$normal <- cbind(sin(comets$i) * sin(comets$om),
                       -sin(comets$i) * cos(comets$om),
                       cos(comets$i))

# Tests to be performed
type_tests <- c("PCvM", "PAD", "PRt")

# Excluding the C/1882 R1-X (Great September comet) records with X = B, C, D
comets_ccf09 <- comets[comets$ccf09, ][-c(13:15), ]

# Sample size
nrow(comets_ccf09)

# Tests for the data in Cuesta-Albertos et al. (2009)
tests_ccf09 <- unif_test(data = comets_ccf09$normal, type = type_tests,
                         p_value = "asymp")
tests_ccf09

---

Craters named by the IUA

Description

Named craters of the Solar System by the Gazetteer of Planetary Nomenclature of the International Astronomical Union (IUA).

Usage

craters

Format

A data frame with 5235 rows and 7 variables:

ID

database ID.

name

name of the crater.

target

name of the celestial body. A factor with 43 levels, such as "Moon", "Venus", or "Europa".

target_type

type of celestial body. A factor with 3 levels: "Planet", "Moon", "Dwarf planet", or "Asteroid".

diameter

diameter of the crater (in km).

theta

longitude angle $\theta \in [0, 2\pi)$ of the crater center.

phi

latitude angle $\phi \in [-\pi/2, \pi/2]$ of the crater center.

Details

"Craters" are understood in the Gazetteer of Planetary Nomenclature as roughly circular depressions resulting from impact or volcanic activity (the geological origin is unspecified).

Be aware that the dataset only contains named craters by the IUA. Therefore, there is likely a high uniform bias on the distribution of craters. Presumably the naming process attempts to cover the planet in a somehow uniform fashion (distant craters are more likely to be named than neighboring craters). Also, there are substantially more craters in the listed bodies than those named by the IUA. See venus and rhea for more detailed and specific crater datasets.

The $(\theta, \phi)$ angles are such their associated planetocentric coordinates are:

$(\cos(\phi) \cos(\theta), \cos(\phi) \sin(\theta), \sin(\phi))',$

with $(0, 0, 1)'$ denoting the north pole.

The script performing the data preprocessing is available at craters.R. The data was retrieved on 2020-05-31.

Source

https://planetarynames.wr.usgs.gov/AdvancedSearch

Examples

# Load data
data("craters")

# Add Cartesian coordinates
craters$X <- cbind(cos(craters$theta) * cos(craters$phi),
                   sin(craters$theta) * cos(craters$phi),
                   sin(craters$phi))

# Tests to be performed
type_tests <- c("PCvM", "PAD", "PRt")

# Tests for Venus and Rhea
unif_test(data = craters$X[craters$target == "Venus", ], type = type_tests,
          p_value = "asymp")
unif_test(data = craters$X[craters$target == "Rhea", ], type = type_tests,
          p_value = "asymp")

---

Distribution and quantile functions from angular function

Description

Numerical computation of the distribution function $F$ and the quantile function $F^{-1}$ for an angular function $f$ in a tangent-normal decomposition. $F^{-1}(x)$ results from the inversion of

$F(x) = \int_{-1}^x \omega_{p - 1}c_f f(z) (1 - z^2)^{(p - 3) / 2} \,\mathrm{d}z$

for $x\in [-1, 1]$, where $c_f$ is a normalizing constant and $\omega_{p - 1}$ is the surface area of $S^{p - 2}$.

Usage

F_from_f(f, p, Gauss = TRUE, N = 320, K = 1000, tol = 1e-06, ...)

F_inv_from_f(f, p, Gauss = TRUE, N = 320, K = 1000, tol = 1e-06, ...)

Arguments

f	angular function defined on $[-1, 1]$. Must be vectorized.
p	integer giving the dimension of the ambient space $R^p$ that contains $S^{p-1}$.
Gauss	use a Gauss–Legendre quadrature rule to integrate $f$ with N nodes? Otherwise, rely on integrate Defaults to TRUE.
N	number of points used in the Gauss–Legendre quadrature. Defaults to 320.
K	number of equispaced points on $[-1, 1]$ used for evaluating $F^{-1}$ and then interpolating. Defaults to 1e3.
tol	tolerance passed to uniroot for the inversion of $F$. Also, passed to integrate's rel.tol and abs.tol if Gauss = FALSE. Defaults to 1e-6.
...	further parameters passed to f.

Details

The normalizing constant $c_f$ is such that $F(1) = 1$. It does not need to be part of f as it is computed internally.

Interpolation is performed by a monotone cubic spline. Gauss = TRUE yields more accurate results, at expenses of a heavier computation.

If f yields negative values, these are silently truncated to zero.

Value

A splinefun object ready to evaluate $F$ or $F^{-1}$, as specified.

Examples

f <- function(x) rep(1, length(x))
plot(F_from_f(f = f, p = 4, Gauss = TRUE), ylab = "F(x)", xlim = c(-1, 1))
plot(F_from_f(f = f, p = 4, Gauss = FALSE), col = 2, add = TRUE,
     xlim = c(-1, 1))
curve(p_proj_unif(x = x, p = 4), col = 3, add = TRUE, n = 300)
plot(F_inv_from_f(f = f, p = 4, Gauss = TRUE), ylab = "F^{-1}(x)")
plot(F_inv_from_f(f = f, p = 4, Gauss = FALSE), col = 2, add = TRUE)
curve(q_proj_unif(u = x, p = 4), col = 3, add = TRUE, n = 300)

---

Gauss–Legendre quadrature

Description

Convenience for computing the nodes $x_k$ and weights $w_k$ of the Gauss–Legendre quadrature formula in $(a, b)$:

$\int_a^b f(x) w(x)\,\mathrm{d}x\approx\sum_{k=1}^N w_k f(x_k).$

.

Usage

Gauss_Legen_nodes(a = -1, b = 1, N = 40L)

Gauss_Legen_weights(a = -1, b = 1, N = 40L)

Arguments

a, b	scalars giving the interval $(a, b)$. Defaults to $(-1, 1)$.
N	number of points used in the Gauss–Legendre quadrature. The following choices are supported: 5, 10, 20, 40, 80, 160, 320, 640, 1280, 2560, and 5120. Defaults to 40.

Details

For $C^\infty$ functions, Gauss–Legendre quadrature can be very efficient. It is exact for polynomials up to degree $2N - 1$.

The nodes and weights up to $N = 80$ were retrieved from NIST and have $10^{-21}$ precision. For $N = 160$ onwards, the nodes and weights were computed with the gauss.quad function from the statmod package (Smyth, 1998), and have $10^{-15}$ precision.

Value

A matrix of size c(N, 1) with the nodes $x_k$ (Gauss_Legen_nodes) or the corresponding weights $w_k$ (Gauss_Legen_weights).

References

NIST Digital Library of Mathematical Functions. Release 1.0.20 of 2018-09-15. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, and B. V. Saunders, eds. https://dlmf.nist.gov/

Smyth, G. K. (1998). Numerical integration. In: Encyclopedia of Biostatistics, P. Armitage and T. Colton (eds.), Wiley, London, pp. 3088-3095.

Examples

## Integration of a smooth function in (1, 10)

# Weights and nodes for integrating
x_k <- Gauss_Legen_nodes(a = 1, b = 10, N = 40)
w_k <- Gauss_Legen_weights(a = 1, b = 10, N = 40)

# Check quadrature
f <- function(x) sin(x) * x^2 - log(x + 1)
integrate(f, lower = 1, upper = 10, rel.tol = 1e-12)
sum(w_k * f(x_k))

# Exact for polynomials up to degree 2 * N - 1
f <- function(x) (((x + 0.5) / 1e3)^5 - ((x - 0.5)/ 5)^4 +
  ((x - 0.25) / 10)^2 + 1)^20
sum(w_k * f(x_k))
integrate(f, lower = -1, upper = 1, rel.tol = 1e-12)

## Integration on (0, pi)

# Weights and nodes for integrating
th_k <- Gauss_Legen_nodes(a = 0, b = pi, N = 40)
w_k <- Gauss_Legen_weights(a = 0, b = pi, N = 40)

# Check quadrature
p <- 4
psi <- function(th) -sin(th / 2)
w <- function(th) sin(th)^(p - 2)
integrate(function(th) psi(th) * w(th), lower = 0, upper = pi,
          rel.tol = 1e-12)
sum(w_k * psi(th_k) * w(th_k))

# Integral with Gegenbauer polynomial
k <- 3
C_k <- function(th) drop(Gegen_polyn(theta = th, k = k, p = p))
integrate(function(th) psi(th) * C_k(th) * w(th), lower = 0, upper = pi,
          rel.tol = 1e-12)
th_k <- drop(Gauss_Legen_nodes(a = 0, b = pi, N = 80))
w_k <- drop(Gauss_Legen_weights(a = 0, b = pi, N = 80))
sum(w_k * psi(th_k) * C_k(th_k) * w(th_k))

---

Gegenbauer polynomials and coefficients

Description

The Gegenbauer polynomials $\{C_k^{(\lambda)}(x)\}_{k = 0}^\infty$ form a family of orthogonal polynomials on the interval $[-1, 1]$ with respect to the weight function $(1 - x^2)^{\lambda - 1/2}$, for $\lambda > -1/2$, $\lambda \neq 0$. They usually appear when dealing with functions defined on $S^{p-1} := \{{\bf x} \in R^p : ||{\bf x}|| = 1\}$ with index $\lambda = p / 2 - 1$.

The Gegenbauer polynomials are somehow simpler to evaluate for $x = \cos(\theta)$, with $\theta \in [0, \pi]$. This simplifies also the connection with the Chebyshev polynomials $\{T_k(x)\}_{k = 0}^\infty$, which admit the explicit expression $T_k(\cos(\theta)) = \cos(k\theta)$. The Chebyshev polynomials appear as the limit of the Gegenbauer polynomials (divided by $\lambda$) when $\lambda$ goes to $0$, so they can be regarded as the extension by continuity of $\{C_k^{(p/2 - 1)}(x)\}_{k = 0}^\infty$ to the case $p = 2$.

For a reasonably smooth function $\psi$ defined on $[0, \pi]$,

$\psi(\theta) = \sum_{k = 0}^\infty b_{k, p} C_k^{(p/2 - 1)}(\cos(\theta)),$

provided that the coefficients

$b_{k, p} := \frac{1}{c_{k, p}} \int_0^\pi \psi(\theta) C_k^{(p/2 - 1)}(\cos(\theta)) (\sin(\theta))^{p - 2}\,\mathrm{d}\theta$

are finite, where the normalizing constants are

$c_{k, p} := \int_0^\pi (C_k^{(p/2 - 1)}(\cos(\theta)))^2 (\sin(\theta))^{p - 2} \,\mathrm{d}\theta.$

The (squared) "Gegenbauer norm" of $\psi$ is

$\|\psi\|_{G, p}^2 := \int_0^\pi \psi(\theta)^2 C_k^{(p/2 - 1)}(\cos(\theta)) (\sin(\theta))^{p - 2}\,\mathrm{d}\theta.$

The previous expansion can be generalized for a 2-dimensional function $\psi$ defined on $[0, \pi] \times [0, \pi]$:

$\psi(\theta_1, \theta_2) = \sum_{k = 0}^\infty \sum_{m = 0}^\infty b_{k, m, p} C_k^{(p/2 - 1)}(\cos(\theta_1)) C_k^{(p/2 - 1)}(\cos(\theta_2)),$

with coefficients

$b_{k, m, p} := \frac{1}{c_{k, p} c_{m, p}} \int_0^\pi\int_0^\pi \psi(\theta_1, \theta_2) C_k^{(p/2 - 1)}(\cos(\theta_1)) C_k^{(p/2 - 1)}(\cos(\theta_2)) (\sin(\theta_1))^{p - 2} (\sin(\theta_2))^{p - 2}\,\mathrm{d}\theta_1\,\mathrm{d}\theta_2.$

The (squared) "Gegenbauer norm" of $\psi$ is

$\|\psi\|_{G, p}^2 := \int_0^\pi\int_0^\pi \psi(\theta_1, \theta_2)^2 C_k^{(p/2 - 1)}(\cos(\theta_1)) C_k^{(p/2 - 1)}(\cos(\theta_2)) (\sin(\theta_1))^{p - 2} (\sin(\theta_2))^{p - 2} \,\mathrm{d}\theta_1\,\mathrm{d}\theta_2.$

Usage

Gegen_polyn(theta, k, p)

Gegen_coefs(k, p, psi, Gauss = TRUE, N = 320, normalize = TRUE,
  only_const = FALSE, tol = 1e-06, ...)

Gegen_series(theta, coefs, k, p, normalize = TRUE)

Gegen_norm(coefs, k, p, normalize = TRUE, cumulative = FALSE)

Gegen_polyn_2d(theta_1, theta_2, k, m, p)

Gegen_coefs_2d(k, m, p, psi, Gauss = TRUE, N = 320, normalize = TRUE,
  only_const = FALSE, tol = 1e-06, ...)

Gegen_series_2d(theta_1, theta_2, coefs, k, m, p, normalize = TRUE)

Gegen_norm_2d(coefs, k, m, p, normalize = TRUE)

Arguments

theta, theta_1, theta_2	vectors with values in $[0, \pi]$.
k, m	vectors with the orders of the Gegenbauer polynomials. Must be integers larger or equal than 0.
p	integer giving the dimension of the ambient space $R^p$ that contains $S^{p-1}$.
psi	function defined in $[0, \pi]$ and whose Gegenbauer coefficients are to be computed. Must be vectorized. For Gegen_coefs_2d, it must return a matrix of size c(length(theta_1), length(theta_2)).
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.
normalize	consider normalized coefficients (divided by $c_{k, p}$)? Defaults to TRUE.
only_const	return only the normalizing constants $c_{k, p}$? Defaults to FALSE.
tol	tolerance passed to integrate's rel.tol and abs.tol if Gauss = FALSE. Defaults to 1e-6.
...	further arguments to be passed to psi.
coefs	for Gegen_series and Gegen_norm, a vector of coefficients $b_{k, p}$ with length length(k). For Gegen_series_2d and Gegen_norm_2d, a matrix of coefficients $b_{k, m, p}$ with size c(length(k), length(m)). The order of the coefficients is given by k and m.
cumulative	return the cumulative norm for increasing truncation of the series? Defaults to FALSE.

Details

The Gegen_polyn function is a wrapper to the functions gegenpoly_n and gegenpoly_array in the gsl-package, which they interface the functions defined in the header file gsl_sf_gegenbauer.h (documented here) of the GNU Scientific Library.

Note that the function Gegen_polyn computes the regular unnormalized Gegenbauer polynomials.

For the case $p = 2$, the Chebyshev polynomials are considered.

Value

• Gegen_polyn: a matrix of size c(length(theta), length(k)) containing the evaluation of the length(k) Gegenbauer polynomials at theta.

• Gegen_coefs: a vector of size length(k) containing the coefficients $b_{k, p}$.

• Gegen_series: the evaluation of the truncated series expansion, a vector of size length(theta).

• Gegen_norm: the Gegenbauer norm of the truncated series, a scalar if cumulative = FALSE, otherwise a vector of size length(k).

• Gegen_polyn_2d: a 4-dimensional array of size c(length(theta_1), length(theta_2), length(k), length(m)) containing the evaluation of the length(k) * length(m) 2-dimensional Gegenbauer polynomials at the bivariate grid spanned by theta_1 and theta_2.

• Gegen_coefs_2d: a matrix of size c(length(k), length(m)) containing the coefficients $b_{k, m, p}$.

• Gegen_series_2d: the evaluation of the truncated series expansion, a matrix of size c(length(theta_1), length(theta_2)).

• Gegen_norm_2d: the 2-dimensional Gegenbauer norm of the truncated series, a scalar.

References

Galassi, M., Davies, J., Theiler, J., Gough, B., Jungman, G., Alken, P., Booth, M., and Rossi, F. (2009) GNU Scientific Library Reference Manual. Network Theory Ltd. http://www.gnu.org/software/gsl/

NIST Digital Library of Mathematical Functions. Release 1.0.20 of 2018-09-15. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller, and B. V. Saunders, eds. https://dlmf.nist.gov/

Examples

## Representation of Gegenbauer polynomials (Chebyshev polynomials for p = 2)

th <- seq(0, pi, l = 500)
k <- 0:3
old_par <- par(mfrow = c(2, 2))
for (p in 2:5) {
  matplot(th, t(Gegen_polyn(theta = th, k = k, p = p)), lty = 1,
          type = "l", main = substitute(p == d, list(d = p)),
          axes = FALSE, xlab = expression(theta), ylab = "")
  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()
  mtext(text = expression({C[k]^{p/2 - 1}}(cos(theta))), side = 2,
        line = 2, cex = 0.75)
  legend("bottomleft", legend = paste("k =", k), lwd = 2, col = seq_along(k))
}
par(old_par)

## Coefficients and series in p = 2

# Function in [0, pi] to be projected in Chebyshev polynomials
psi <- function(th) -sin(th / 2)

# Coefficients
p <- 2
k <- 0:4
(coefs <- Gegen_coefs(k = k, p = p, psi = psi))

# Series
plot(th, psi(th), type = "l", axes = FALSE, xlab = expression(theta),
      ylab = "", ylim = c(-1.25, 0))
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()
col <- viridisLite::viridis(length(coefs))
for (i in seq_along(coefs)) {
  lines(th, Gegen_series(theta = th, coefs = coefs[1:(i + 1)], k = 0:i,
                         p = p), col = col[i])
}
lines(th, psi(th), lwd = 2)

## Coefficients and series in p = 3

# Function in [0, pi] to be projected in Gegenbauer polynomials
psi <- function(th) tan(th / 3)

# Coefficients
p <- 3
k <- 0:10
(coefs <- Gegen_coefs(k = k, p = p, psi = psi))

# Series
plot(th, psi(th), type = "l", axes = FALSE, xlab = expression(theta),
      ylab = "", ylim = c(0, 2))
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()
col <- viridisLite::viridis(length(coefs))
for (i in seq_along(coefs)) {
  lines(th, Gegen_series(theta = th, coefs = coefs[1:(i + 1)], k = 0:i,
                         p = p), col = col[i])
}
lines(th, psi(th), lwd = 2)

## Surface representation

# Surface in [0, pi]^2 to be projected in Gegenbauer polynomials
p <- 3
psi <- function(th_1, th_2) A_theta_x(theta = th_1, x = cos(th_2),
                                      p = p, as_matrix = TRUE)

# Coefficients
k <- 0:20
m <- 0:10
coefs <- Gegen_coefs_2d(k = k, m = m, p = p, psi = psi)

# Series
th <- seq(0, pi, l = 100)
col <- viridisLite::viridis(20)
old_par <- par(mfrow = c(2, 2))
image(th, th, A_theta_x(theta = th, x = cos(th), p = p), axes = FALSE,
      col = col, zlim = c(0, 1), xlab = expression(theta[1]),
      ylab = expression(theta[2]), main = "Original")
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, at = c(0, pi / 4, pi / 2, 3 * pi / 4, pi),
     labels = expression(0, pi / 4, pi / 2, 3 * pi / 4, pi))
box()
for(K in c(5, 10, 20)) {
  A <- Gegen_series_2d(theta_1 = th, theta_2 = th,
                       coefs = coefs[1:(K + 1), ], k = 0:K, m = m, p = p)
  image(th, th, A, axes = FALSE, col = col, zlim = c(0, 1),
        xlab = expression(theta[1]), ylab = expression(theta[2]),
        main = paste(K, "x", m[length(m)], "coefficients"))
  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, at = c(0, pi / 4, pi / 2, 3 * pi / 4, pi),
       labels = expression(0, pi / 4, pi / 2, 3 * pi / 4, pi))
  box()
}
par(old_par)

---

(Hyper)spherical harmonics

Description

Computation of a certain explicit representation of (hyper)spherical harmonics on $S^{p-1}:=\{{\bf x}\in R^p:||{\bf x}||=1\}$, $p\ge 2$. Details are available in García-Portugués et al. (2024).

Usage

g_i_k(x, i = 1, k = 1, m = NULL, show_m = FALSE)

Arguments

x	locations in $S^{p-1}$ to evaluate $g_{i,k}$. Either a matrix of size c(nx, p) or a vector of size p. Normalized internally if required (with a warning message).
i, k	alternative indexing to refer to the i-th (hyper)spherical harmonic of order k. i is a positive integer smaller than d_p_k and k is a non-negative integer.
m	(hyper)spherical harmonic index, as used in Proposition 3.1. The index is computed internally from i and k. Defaults to NULL.
show_m	flag to print m if computed internally when m = NULL.

Details

The implementation uses Proposition 3.1 in García-Portugués et al. (2024), which adapts Theorem 1.5.1 in Dai and Xu (2013) with the correction of typos in the normalizing constant $h_\alpha$ and in the definition of the function $g_\alpha$ of the latter theorem.

Value

A vector of size nrow(x).

References

Dai, F. and Xu, Y. (2013). Approximation Theory and Harmonic Analysis on Spheres and Balls. Springer, New York. doi:10.1007/978-1-4614-6660-4

García-Portugués, E., Paindaveine, D., and Verdebout, T. (2024). On a class of Sobolev tests for symmetry of directions, their detection thresholds, and asymptotic powers. arXiv:2108.09874v2. doi:10.48550/arXiv.2108.09874

Examples

n <- 3e3
old_par <- par(mfrow = c(2, 3))
k <- 2
for (i in 1:d_p_k(p = 3, k = k)) {
  X <- r_unif_sph(n = n, p = 3, M = 1)[, , 1]
  col <- rainbow(n)[rank(g_i_k(x = X, k = k, i = i, show_m = TRUE))]
  scatterplot3d::scatterplot3d(X[, 1], X[, 2], X[, 3], color = col,
                               axis = FALSE, pch = 19)
}
for (k in 0:5) {
  X <- r_unif_sph(n = n, p = 3, M = 1)[, , 1]
  col <- rainbow(n)[rank(g_i_k(x = X, k = k, i = 1, show_m = TRUE))]
  scatterplot3d::scatterplot3d(X[, 1], X[, 2], X[, 3], color = col,
                               axis = FALSE, pch = 19)
}
par(old_par)

---

Monte Carlo integration of functions on the (hyper)sphere

Description

Monte Carlo approximation of the integral

$\int_{S^{p-1}}f(x)\,\mathrm{d}x$

of a function $f:S^{p-1} \rightarrow R$ defined on the (hyper)sphere $S^{p-1}:=\{{\bf x}\in R^p:||{\bf x}||=1\}$, $p\ge 2$.

Usage

int_sph_MC(f, p, M = 10000, cores = 1, chunks = ceiling(M/1000),
  seeds = NULL, ...)

Arguments

f	function to be integrated. Its first argument must be the (hyper)sphere position. Must be vectorized and return a vector of size nrow(x) for a matrix input x. See examples.
p	integer giving the dimension of the ambient space $R^p$ that contains $S^{p-1}$.
M	number of Monte Carlo samples. Defaults to 1e4.
cores	number of cores to perform the integration. Defaults to 1.
chunks	number of chunks to split the M Monte Carlo samples. Useful for parallelizing the integration in chunks tasks containing ceiling(M / chunks) replications. Useful also for avoiding memory bottlenecks when M is large. Defaults to ceiling(M / 1e3).
seeds	if provided, a vector of size chunks for fixing the seeds on each of the simulation chunks (useful for reproducing parallel simulations). Specifically, for k in 1:chunks, seeds are set as set.seed(seeds[k], kind = "Mersenne-Twister") in each chunk. Defaults to NULL (no seed setting is done).
...	optional arguments to be passed to f or to foreach (for example, .export to export global variables or other functions to the foreach environment).

Details

It is possible to have a progress bar if int_sph_MC is wrapped with progressr::with_progress or if progressr::handlers(global = TRUE) is invoked (once) by the user. See the examples below. The progress bar is updated with the number of finished chunks.

Value

A scalar with the approximate integral.

Examples

## Sequential simulation

# Vectorized functions to be integrated
x1 <- function(x) x[, 1]
quad <- function(x, a = 0) a + rowSums(x^4)

# Approximate \int_{S^{p-1}} x_1 dx = 0
int_sph_MC(f = x1, p = 3, M = 1e4, chunks = 2)

# Approximate \int_{S^{p-1}} (a + \sum_i x_i^4) dx
int_sph_MC(f = quad, p = 2, M = 1e4, a = 0, chunks = 2)

# Compare with Gauss--Legendre integration on S^2
th_k <- Gauss_Legen_nodes(a = 0, b = 2 * pi, N = 40)
w_k <- Gauss_Legen_weights(a = 0, b = 2 * pi, N = 40)
sum(w_k * quad(cbind(cos(th_k), sin(th_k)), a = 1))

## Parallel simulation with a progress bar

# Define a progress bar
require(progress)
require(progressr)
handlers(handler_progress(
  format = paste("(:spin) [:bar] :percent Iter: :current/:total Rate:",
                 ":tick_rate iter/sec ETA: :eta Elapsed: :elapsedfull"),
  clear = FALSE))
# Call int_sph_MC() within with_progress()
with_progress(int_sph_MC(f = x1, p = 3, cores = 2, M = 1e5, chunks = 100))

# Instead of using with_progress() each time, it is more practical to run
# handlers(global = TRUE)
# once to activate progress bars in your R session

---

Local projected alternatives to uniformity

Description

Density and random generation for local projected alternatives to uniformity with densities

$f_{\kappa, \boldsymbol{\mu}}({\bf x}): = \frac{1 - \kappa}{\omega_p} + \kappa f({\bf x}'\boldsymbol{\mu})$

where

$f(z) = \frac{1}{\omega_p}\left\{1 + \sum_{k = 1}^\infty u_{k, p} C_k^{p / 2 - 1}(z)\right\}$

is the angular function controlling the local alternative in a Gegenbauer series, $0\le \kappa \le 1$, $\boldsymbol{\mu}$ is a direction on $S^{p - 1}$, and $\omega_p$ is the surface area of $S^{p - 1}$. The sequence $\{u_{k, p}\}$ is typically such that $u_{k, p} = \left(1 + \frac{2k}{p - 2}\right) b_{k, p}$ for the Gegenbauer coefficients $\{b_{k, p}\}$ of the kernel function of a Sobolev statistic (see the transformation between the coefficients $u_{k, p}$ and $b_{k, p}$).

Also, automatic truncation of the series $\sum_{k = 1}^\infty u_{k, p} C_k^{p / 2 - 1}(z)$ according to the proportion of "Gegenbauer norm" explained.

Usage

f_locdev(z, p, uk)

con_f(f, p, N = 320)

d_locdev(x, mu, f, kappa)

r_locdev(n, mu, f, kappa, F_inv = NULL, ...)

cutoff_locdev(p, K_max = 10000, thre = 0.001, type, Rothman_t = 1/3,
  Pycke_q = 0.5, verbose = FALSE, Gauss = TRUE, N = 320, tol = 1e-06)

Arguments

z	projected evaluation points for $f$, a vector with entries on $[-1, 1]$.
p	integer giving the dimension of the ambient space $R^p$ that contains $S^{p-1}$.
uk	coefficients $u_{k, p}$ associated to the indexes 1:length(uk), a vector.
f	angular function defined on $[-1, 1]$. Must be vectorized.
N	number of points used in the Gauss–Legendre quadrature for computing the Gegenbauer coefficients. Defaults to 320.
x	locations in $S^{p-1}$ to evaluate the density. Either a matrix of size c(nx, p) or a vector of length p. Normalized internally if required (with a warning message).
mu	a unit norm vector of size p giving the axis of rotational symmetry.
kappa	the strength of the local alternative, between 0 and 1.
n	sample size, a positive integer.
F_inv	quantile function associated to $f$. Computed by F_inv_from_f if NULL (default).
...	further parameters passed to F_inv_from_f.
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.
type	name of the Sobolev statistic, using the naming from avail_cir_tests and avail_sph_tests.
Rothman_t	$t$ parameter for the Rothman test, a real in $(0, 1)$. Defaults to 1 / 3.
Pycke_q	$q$ parameter for the Pycke "$q$-test", a real in $(0, 1)$. Defaults to 1 / 2.
verbose	output information about the truncation (TRUE or 1) and a diagnostic plot (2)? Defaults to FALSE.
Gauss	use a Gauss–Legendre quadrature rule of N nodes in the computation of the Gegenbauer coefficients? Otherwise, call integrate. Defaults to TRUE.
tol	tolerance passed to integrate's rel.tol and abs.tol if Gauss = FALSE. Defaults to 1e-6.

Details

See the definitions of local alternatives in Prentice (1978) and in García-Portugués et al. (2023).

The truncation of $\sum_{k = 1}^\infty u_{k, p} C_k^{p / 2 - 1}(z)$ is done to the first K_max terms and then up to the index such that the first terms leave unexplained the proportion thre of the norm of the whole series. Setting thre = 0 truncates to K_max terms exactly. If the series only contains odd or even non-zero terms, then only K_max / 2 addends are effectively taken into account in the first truncation.

Value

• f_locdev: angular function evaluated at x, a vector.

• con_f: normalizing constant $c_f$ of $f$, a scalar.

• d_locdev: density function evaluated at x, a vector.

• r_locdev: a matrix of size c(n, p) containing a random sample from the density $f_{\kappa, \boldsymbol{\mu}}$.

• cutoff_locdev: vector of coefficients $\{u_{k, p}\}$ automatically truncated according to K_max and thre (see 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. doi:10.3150/21-BEJ1454.

Prentice, M. J. (1978). On invariant tests of uniformity for directions and orientations. The Annals of Statistics, 6(1):169–176. doi:10.1214/aos/1176344075

Examples

## Local alternatives diagnostics

loc_alt_diagnostic  <- function(p, type, thre = 1e-3, K_max = 1e3) {

  # Coefficients of the alternative
  uk <- cutoff_locdev(K_max = K_max, p = p, type = type, thre = thre,
                      N = 640)

  old_par <- par(mfrow = c(2, 2))

  # Construction of f
  z <- seq(-1, 1, l = 1e3)
  f <- function(z) f_locdev(z = z, p = p, uk = uk)
  plot(z, f(z), type = "l", xlab = expression(z), ylab = expression(f(z)),
       main = paste0("Local alternative f, ", type, ", p = ", p), log = "y")

  # Projected density on [-1, 1]
  f_proj <- function(z) rotasym::w_p(p = p - 1) * f(z) *
    (1 - z^2)^((p - 3) / 2)
  plot(z, f_proj(z), type = "l", xlab = expression(z),
       ylab = expression(omega[p - 1] * f(z) * (1 - z^2)^{(p - 3) / 2}),
       main = paste0("Projected density, ", type, ", p = ", p), log = "y",
       sub = paste("Integral:", round(con_f(f = f, p = p), 4)))

  # Quantile function for projected density
  mu <- c(rep(0, p - 1), 1)
  F_inv <- F_inv_from_f(f = f, p = p, K = 5e2)
  plot(F_inv, xlab = expression(x), ylab = expression(F^{-1}*(x)),
       main = paste0("Quantile function, ", type, ", p = ", p))

  # Sample from the alternative and plot the projected sample
  n <- 5e4
  samp <- r_locdev(n = n, mu = mu, f = f, kappa = 1, F_inv = F_inv)
  plot(z, f_proj(z), col = 2, type = "l",
       main = paste0("Simulated projected data, ", type, ", p = ", p),
       ylim = c(0, 1.75))
  hist(samp %*% mu, freq = FALSE, breaks = seq(-1, 1, l = 50), add = TRUE)

  par(old_par)

}

## Local alternatives for the PCvM test

loc_alt_diagnostic(p = 2, type = "PCvM")
loc_alt_diagnostic(p = 3, type = "PCvM")
loc_alt_diagnostic(p = 4, type = "PCvM")
loc_alt_diagnostic(p = 5, type = "PCvM")
loc_alt_diagnostic(p = 11, type = "PCvM")

## Local alternatives for the PAD test

loc_alt_diagnostic(p = 2, type = "PAD")
loc_alt_diagnostic(p = 3, type = "PAD")
loc_alt_diagnostic(p = 4, type = "PAD")
loc_alt_diagnostic(p = 5, type = "PAD")
loc_alt_diagnostic(p = 11, type = "PAD")

## Local alternatives for the PRt test

loc_alt_diagnostic(p = 2, type = "PRt")
loc_alt_diagnostic(p = 3, type = "PRt")
loc_alt_diagnostic(p = 4, type = "PRt")
loc_alt_diagnostic(p = 5, type = "PRt")
loc_alt_diagnostic(p = 11, type = "PRt")

---