r_alt function

Sample non-uniformly distributed spherical data

Simple simulation of prespecified non-uniform spherical distributions: von Mises--Fisher (vMF), Mixture of vMF (MvMF), Angular Central Gaussian (ACG), Small Circle (SC), Watson (W), Cauchy-like (C), Mixture of Cauchy-like (MC), or Uniform distribution with Antipodal-Dependent observations (UAD).

r_alt(n, p, M = 1, alt = "vMF", mu = c(rep(0, p - 1), 1), kappa = 1,
  nu = 0.5, F_inv = NULL, K = 1000, axial_mix = TRUE)

Arguments

• n: sample size.
• p: integer giving the dimension of the ambient space $R^p$ that contains $S^{p-1}$.
• M: number of samples of size n. Defaults to 1.
• alt: alternative, must be "vMF", "MvMF", "ACG", "SC", "W", "C", "MC", or "UAD". See details below.
• mu: location parameter for "vMF", "SC", "W", and "C". Defaults to c(rep(0, p - 1), 1).
• kappa: non-negative parameter measuring the strength of the deviation with respect to uniformity (obtained with $\kappa = 0$).
• nu: projection along $e_p$ controlling the modal strip of the small circle distribution. Must be in (-1, 1). Defaults to 0.5.
• F_inv: quantile function returned by F_inv_from_f. Used for "SC", "W", and "C". Computed by internally if NULL (default).
• K: number of equispaced points on $[-1, 1]$ used for evaluating $F^{-1}$ and then interpolating. Defaults to 1e3.
• axial_mix: use a mixture of von Mises--Fisher or Cauchy-like that is axial (i.e., symmetrically distributed about the origin)? Defaults to TRUE.

Returns

An array of size c(n, p, M) with M random samples of size n of non-uniformly-generated directions on $S^{p-1}$.

Details

The parameter kappa is used as $\kappa$ in the following distributions:

• "vMF": von Mises--Fisher distribution with concentration $\kappa$ and directional mean $\boldsymbol{\mu}$.
• "MvMF": equally-weighted mixture of $p$ von Mises--Fisher distributions with common concentration $\kappa$ and directional means $±e_1, \ldots, ±e_p$ if axial_mix = TRUE. If axial_mix = FALSE, then only means with positive signs are considered.
• "ACG": Angular Central Gaussian distribution with diagonal shape matrix with diagonal given by

• "SC": Small Circle distribution with axis mean $\boldsymbol{\mu}$ and concentration $\kappa$ about the projection along the mean, $\nu$.

• "W": Watson distribution with axis mean $\boldsymbol{\mu}$

  and concentration $\kappa$. The Watson distribution is a particular case of the Bingham distribution.

• "C": Cauchy-like distribution with directional mode $\boldsymbol{\mu}$ and concentration $\kappa = \rho / (1 - \rho^2)$. The circular Wrapped Cauchy distribution is a particular case of this Cauchy-like distribution.

• "MC": equally-weighted mixture of $p$ Cauchy-like distributions with common concentration $\kappa$ and directional means $±e_1, \ldots, ±e_p$ if axial_mix = TRUE. If axial_mix = FALSE, then only means with positive signs are considered.

The alternative "UAD" generates a sample formed by $\lceil n/2\rceil$ observations drawn uniformly on $S^{p-1}$

and the remaining observations drawn from a uniform spherical cap distribution of angle $\pi-\kappa$ about each of the $\lceil n/2\rceil$ observations (see unif_cap). Hence, kappa = 0 corresponds to a spherical cap covering the whole sphere and kappa = pi is a one-point degenerate spherical cap.

Much faster sampling for "SC", "W", "C", and "MC"

is achieved providing F_inv; see examples.

Examples

## Simulation with p = 2

p <- 2
n <- 50
kappa <- 20
nu <- 0.5
angle <- pi / 10
rho <- ((2 * kappa + 1) - sqrt(4 * kappa + 1)) / (2 * kappa)
F_inv_SC_2 <- F_inv_from_f(f = function(z) exp(-kappa * (z - nu)^2), p = 2)
F_inv_W_2 <- F_inv_from_f(f = function(z) exp(kappa * z^2), p = 2)
F_inv_C_2 <- F_inv_from_f(f = function(z) (1 - rho^2) /
                            (1 + rho^2 - 2 * rho * z)^(p / 2), p = 2)
x1 <- r_alt(n = n, p = p, alt = "vMF", kappa = kappa)[, , 1]
x2 <- r_alt(n = n, p = p, alt = "C", F_inv = F_inv_C_2)[, , 1]
x3 <- r_alt(n = n, p = p, alt = "SC", F_inv = F_inv_SC_2)[, , 1]
x4 <- r_alt(n = n, p = p, alt = "ACG", kappa = kappa)[, , 1]
x5 <- r_alt(n = n, p = p, alt = "W", F_inv = F_inv_W_2)[, , 1]
x6 <- r_alt(n = n, p = p, alt = "MvMF", kappa = kappa)[, , 1]
x7 <- r_alt(n = n, p = p, alt = "MC", kappa = kappa)[, , 1]
x8 <- r_alt(n = n, p = p, alt = "UAD", kappa = 1 - angle)[, , 1]
r <- runif(n, 0.95, 1.05) # Radius perturbation to improve visualization
plot(r * x1, pch = 16, xlim = c(-1.1, 1.1), ylim = c(-1.1, 1.1), col = 1)
points(r * x2, pch = 16, col = 2)
points(r * x3, pch = 16, col = 3)
plot(r * x4, pch = 16, xlim = c(-1.1, 1.1), ylim = c(-1.1, 1.1), col = 1)
points(r * x5, pch = 16, col = 2)
points(r * x6, pch = 16, col = 3)
points(r * x7, pch = 16, col = 4)
col <- rep(rainbow(n / 2), 2)
plot(r * x8, pch = 16, xlim = c(-1.1, 1.1), ylim = c(-1.1, 1.1), col = col)
for (i in seq(1, n, by = 2)) lines((r * x8)[i + 0:1, ], col = col[i])

## Simulation with p = 3

n <- 50
p <- 3
kappa <- 20
angle <- pi / 10
nu <- 0.5
rho <- ((2 * kappa + 1) - sqrt(4 * kappa + 1)) / (2 * kappa)
F_inv_SC_3 <- F_inv_from_f(f = function(z) exp(-kappa * (z - nu)^2), p = 3)
F_inv_W_3 <- F_inv_from_f(f = function(z) exp(kappa * z^2), p = 3)
F_inv_C_3 <- F_inv_from_f(f = function(z) (1 - rho^2) /
                            (1 + rho^2 - 2 * rho * z)^(p / 2), p = 3)
x1 <- r_alt(n = n, p = p, alt = "vMF", kappa = kappa)[, , 1]
x2 <- r_alt(n = n, p = p, alt = "C", F_inv = F_inv_C_3)[, , 1]
x3 <- r_alt(n = n, p = p, alt = "SC", F_inv = F_inv_SC_3)[, , 1]
x4 <- r_alt(n = n, p = p, alt = "ACG", kappa = kappa)[, , 1]
x5 <- r_alt(n = n, p = p, alt = "W", F_inv = F_inv_W_3)[, , 1]
x6 <- r_alt(n = n, p = p, alt = "MvMF", kappa = kappa)[, , 1]
x7 <- r_alt(n = n, p = p, alt = "MC", kappa = kappa)[, , 1]
x8 <- r_alt(n = n, p = p, alt = "UAD", kappa = 1 - angle)[, , 1]
s3d <- scatterplot3d::scatterplot3d(x1, pch = 16, xlim = c(-1.1, 1.1),
                                    ylim = c(-1.1, 1.1), zlim = c(-1.1, 1.1))
s3d$points3d(x2, pch = 16, col = 2)
s3d$points3d(x3, pch = 16, col = 3)
s3d <- scatterplot3d::scatterplot3d(x4, pch = 16, xlim = c(-1.1, 1.1),
                                    ylim = c(-1.1, 1.1), zlim = c(-1.1, 1.1))
s3d$points3d(x5, pch = 16, col = 2)
s3d$points3d(x6, pch = 16, col = 3)
s3d$points3d(x7, pch = 16, col = 4)
col <- rep(rainbow(n / 2), 2)
s3d <- scatterplot3d::scatterplot3d(x8, pch = 16, xlim = c(-1.1, 1.1),
                                    ylim = c(-1.1, 1.1), zlim = c(-1.1, 1.1),
                                    color = col)
for (i in seq(1, n, by = 2)) s3d$points3d(x8[i + 0:1, ], col = col[i],
                                          type = "l")