sphunif (version 1.4.3)

r_alt: Sample non-uniformly distributed spherical data

Description

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

Usage

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)

Value

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

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 \({\bf 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.

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 \(\pm{\bf e}_1, \ldots, \pm{\bf 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 $$(1, \ldots, 1, 1 + \kappa) / (p + \kappa).$$

• "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 \(\pm{\bf e}_1, \ldots, \pm{\bf 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")