Random sampling from the Poisson kernel-based density

Giovanni Saraceno

library(QuadratiK)

In this example, we generate observations from the Poisson kernel-based distribution on the sphere, \(S^{d-1}\). We consider mean direction \(\mu=(0,0,1)\), \(d=3\) and the concentration parameter is \(\rho = 0.8\).

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

We sampled \(n=1000\) observations for each method available.

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 function returns a list with the following components

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()
## wgl 
##   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()