r_kern_polysph: Sample kernel-distributed polyspherical data

In polykde: Polyspherical Kernel Density Estimation

Sample kernel-distributed polyspherical data

Description

Simulates from the distribution defined by a kernel on the polysphere \mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}.

Usage

r_kern_polysph(n, d, mu, h, kernel = 1, kernel_type = 1, k = 10,
  intrinsic = FALSE, norm_mu = FALSE)

Arguments

n	sample size.
d	vector of size r with dimensions.
mu	a vector of size sum(d) + r with the concatenated means that define the center of the kernel.
h	vector of size r with bandwidths.
kernel	kernel employed: 1 for von Mises–Fisher (default); 2 for Epanechnikov; 3 for softplus.
kernel_type	type of kernel employed: 1 for product kernel (default); 2 for spherically symmetric kernel.
k	softplus kernel parameter. Defaults to 10.0.
intrinsic	use the intrinsic distance, instead of the extrinsic-chordal distance, in the kernel? Defaults to FALSE.
norm_mu	ensure a normalization of mu? Defaults to FALSE.

Details

Simulation for non-von Mises–Fisher spherically symmetric kernels is done by acceptance-rejection from a von Mises–Fisher proposal distribution.

Value

A matrix of size c(n, sum(d) + r) with the sample.

Examples

# Simulate kernels in (S^1)^2
n <- 1e3
h <- c(1, 1)
d <- c(1, 1)
mu <- rep(DirStats::to_cir(pi), 2)
samp_ker <- function(kernel, kernel_type, col, main) {
  data <- r_kern_polysph(n = n, d = d, mu = mu, h = h, kernel = kernel,
                         kernel_type = kernel_type)
  ang <- cbind(DirStats::to_rad(data[, 1:2]),
               DirStats::to_rad(data[, 3:4]))
  plot(ang, xlim = c(0, 2 * pi), ylim = c(0, 2 * pi), pch = 16, cex = 0.25,
       col = col, xlab = expression(Theta[1]), ylab = expression(Theta[2]),
       main = main)
}
old_par <- par(mfcol = c(2, 3))
samp_ker(kernel = 2, kernel_type = 2, col = 1, main = "Epa sph. symmetric")
samp_ker(kernel = 2, kernel_type = 1, col = 2, main = "Epa product")
samp_ker(kernel = 3, kernel_type = 2, col = 1, main = "Sfp sph. symmetric")
samp_ker(kernel = 3, kernel_type = 1, col = 2, main = "Sfp product")
samp_ker(kernel = 1, kernel_type = 2, col = 1, main = "vMF sph. symmetric")
samp_ker(kernel = 1, kernel_type = 1, col = 2, main = "vMF product")
par(old_par)

# Simulate kernels in (S^2)^2
n <- 1e3
h <- c(0.2, 0.6)
d <- c(2, 2)
mu <- c(c(0, 0, 1), c(0, -1, 0))
samp_ker <- function(kernel, kernel_type, main) {
  data <- r_kern_polysph(n = n, d = d, mu = mu, h = h, kernel = kernel,
                         kernel_type = kernel_type)
  scatterplot3d::scatterplot3d(rbind(data[, 1:3], data[, 4:6]),
                               xlim = c(-1, 1), ylim = c(-1, 1),
                               zlim = c(-1, 1), pch = 16, xlab = "",
                               ylab = "", zlab = "", cex.symbols = 0.5,
       color = rep(viridis::viridis(n)[rank(data[, 3])], 2), main = main)
}
old_par <- par(mfcol = c(2, 3))
samp_ker(kernel = 2, kernel_type = 2, main = "Epa sph. symmetric")
samp_ker(kernel = 2, kernel_type = 1, main = "Epa product")
samp_ker(kernel = 3, kernel_type = 2, main = "Sfp sph. symmetric")
samp_ker(kernel = 3, kernel_type = 1, main = "Sfp product")
samp_ker(kernel = 1, kernel_type = 2, main = "vMF sph. symmetric")
samp_ker(kernel = 1, kernel_type = 1, main = "vMF product")
par(old_par)

# Plot simulated data
n <- 1e3
h <- c(1, 1)
d <- c(2, 2)
samp_ker <- function(kernel, kernel_type, col, main) {
  X <- r_kern_polysph(n = n, d = d, mu = mu, h = h, kernel = kernel,
                      kernel_type = kernel_type)
  S <- cbind((1 - X[, 1:3] %*% mu[1:3]) / h[1]^2,
             (1 - X[, 4:6] %*% mu[4:6]) / h[2]^2)
  plot(S, xlim = c(0, 2 / h[1]^2), ylim = c(0, 2 / h[2]^2), pch = 16,
       cex = 0.25, col = col, xlab = expression(t[1]),
       ylab = expression(t[2]), main = main)
  t_grid <- seq(0, 2 / min(h)^2, l = 100)
  gr <- as.matrix(expand.grid(t_grid, t_grid))
  if (kernel_type == "1") {

    dens <- prod(c_kern(h = h, d = d, kernel = kernel, kernel_type = 1)) *
      L(gr[, 1], kernel = kernel) * L(gr[, 2], kernel = kernel)

  } else if (kernel_type == "2") {

    dens <- c_kern(h = h, d = d, kernel = kernel, kernel_type = 2) *
      L(gr[, 1] + gr[, 2], kernel = kernel)

  }
  dens <- matrix(dens, nrow = length(t_grid), ncol = length(t_grid))
  contour(t_grid, t_grid, dens, add = TRUE, col = col,
          levels = seq(0, 0.2, l = 41))

}
old_par <- par(mfcol = c(2, 3))
samp_ker(kernel = 2, kernel_type = 2, col = 1, main = "Epa sph. symmetric")
samp_ker(kernel = 2, kernel_type = 1, col = 2, main = "Epa product")
samp_ker(kernel = 3, kernel_type = 2, col = 1, main = "Sfp sph. symmetric")
samp_ker(kernel = 3, kernel_type = 1, col = 2, main = "Sfp product")
samp_ker(kernel = 1, kernel_type = 2, col = 1, main = "vMF sph. symmetric")
samp_ker(kernel = 1, kernel_type = 1, col = 2, main = "vMF product")
par(old_par)

polykde documentation built on April 16, 2025, 1:11 a.m.