Random sampling from the Poisson kernel-based density

library(QuadratiK)

In this example, we generate observations from the Poisson kernel-based distribution on the sphere, S^d − 1. Given a vector μ ∈ 𝒮^d − 1, where 𝒮^d − 1 = {x ∈ ℝ^d : ||x|| = 1}, and a parameter ρ such that 0 < ρ < 1, the probability density function of a d-variate Poisson kernel-based density is defined by:

$$ f(\mathbf{x}|\rho, \mathbf{\mu}) = \frac{1-\rho^2}{\omega_d ||\mathbf{x} - \rho \mathbf{\mu}||^d}, $$

where μ is a vector orienting the center of the distribution, ρ is a parameter to control the concentration of the distribution around the vector μ, and it is related to the variance of the distribution. Furthermore, ω_d = 2π^d/2[Γ(d/2)]⁻¹ is the surface area of the unit sphere in ℝ^d (see Golzy and Markatou, 2020). When ρ → 0, the Poisson kernel-based density tends to the uniform density on the sphere. Connection of the PKBDs to other distributions are discussed in detail in Golzy and Markatou (2020).

We consider mean direction μ = (0, 0, 1), d = 3 and the concentration parameter is ρ = 0.8.

mu <- c(0,0,1)
d <- 3
n <- 1000
rho <- 0.8

We sampled n = 1000 observations for each method available. Golzy and Markatou (2020) proposed an acceptance-rejection method for simulating data from a PKBD using von Mises-Fisher envelopes (rejvmf method). Furthermore, Sablica, Hornik, and Leydold (2023) proposed new ways for simulating from the PKBD, using angular central Gaussian envelopes (rejacg) or using the projected Saw distributions (rejpsaw).

n <- 1000
set.seed(2468)
# Generate observations using the rejection algorithm with von-Mises 
# distribution envelopes
dat1 <- rpkb(n = n, rho=rho, mu=mu, method="rejvmf")
# Generate observations using the rejection algorithm with angular central 
# Gaussian distribution envelopes
dat2 <- rpkb(n = n, rho=rho, mu=mu, method="rejacg")
# Generate observations using the projected Saw distribution
dat3 <- rpkb(n = n, rho=rho, mu=mu, method="rejpsaw")

The number of observations generated is determined by n for rpkb(). This function returns a list with the matrix of generated observations x, the number of tries numTries, and the number of acceptances numAccepted.

summary(dat1)

##             Length Class  Mode   
## x           3000   -none- numeric
## numAccepted    1   -none- numeric
## numTries       1   -none- numeric

We extract the generated data sets and create a unique data set.

x <- rbind(dat1$x, dat2$x, dat3$x)

Finally, the observations generated through the different methods are compared graphically, by displaying the data points on the sphere colored with respect to the used method.

library(rgl)
# Legend information
classes <- c("rejvmf", "rejacg", "rejpsaw")
# Fix a color for each method
colors <- c("red", "blue", "green")
labels <- factor(rep(colors, each = 1000))
# Element needed for the Legend position in the following plot
offset <- 0.25
# Coordinates for legend placement
legend_x <- max(x[,1]) + offset  
legend_y <- max(x[,2]) + offset
legend_z <- seq(min(x[,3]), length.out = length(classes), by = offset)

open3d()

## null 
##    1

# Create the legend
for (i in seq_along(classes)) {
   text3d(legend_x, legend_y, legend_z[i], texts = classes[i], adj = c(0, 0.5))
   points3d(legend_x-0.1, legend_y, legend_z[i], col = colors[i], size = 5)
}
title3d("", line = 3, cex = 1.5, font=2, add=TRUE)
# Plot the sampled observations colored with respect to the used method
plot3d(x[,1], x[,2], x[,3], col = labels, size = 5, add=TRUE)
title3d("", line = 3, cex = 1.5, font=2, add=TRUE)
# Add a Sphere as background
rgl.spheres(0 , col = "transparent", alpha = 0.2)
# Rotate and zoom the generated 3 dimensional plot to facilitate visualization
view3d(theta = 10, phi = -25, zoom = 0.5)
# rglwidget is added in order to display the generated figure into the html 
# replication file.
rglwidget()

close3d()