# Install the package
install.packages("DirStats")

# Load package
library(DirStats)

# Install the package
library(devtools)
install_github("egarpor/DirStats")

# Load package
library(DirStats)

# Sample from a von Mises--Fisher on S^2
q <- 2
n <- 300
set.seed(42)
samp <- rbind(rotasym::r_vMF(n = n / 3, mu = c(rep(0, q), 1), kappa = 5),
              rotasym::r_vMF(n = n / 3, mu = c(rep(0, q), -1), kappa = 5),
              rotasym::r_vMF(n = n / 3, mu = c(1, rep(0, q)), kappa = 5))

# LCV bandwidth
(h_LCV <- bw_dir_lcv(data = samp)$h_opt)
#> [1] 0.2567438

# LSCV bandwidth
(h_LSCV <- bw_dir_lscv(data = samp)$h_opt)
#> [1] 0.232965

# ROT bandwidth
(h_ROT <- bw_dir_rot(data = samp))
#> [1] 0.4053648

# Mixture fit, for AMI and EMI bandwidth selectors
fit_mix <- bic_vmf_mix(data = samp)

# AMI bandwidth
(h_AMI <- bw_dir_ami(data = samp, fit_mix = fit_mix))
#> [1] 0.2054242

# EMI bandwidth
(h_EMI <- bw_dir_emi(data = samp, fit_mix = fit_mix)$h_opt)
#> [1] 0.22527

# Compute kde
l <- 200
th <- seq(0, 2 * pi, l = l)
phi <- seq(0, pi, l = l / 2)
ang <- expand.grid(th = th, phi = phi)
x <- to_sph(th = ang$th, ph = ang$phi)
kde <- kde_dir(x = x, data = samp, h = h_EMI)

# Visualization
contour(x = th, y = phi, z = matrix(kde, nrow = l, ncol = l / 2),
        col = viridisLite::viridis(15), 
        levels = round(seq(min(kde), max(kde), length.out = 15), 4), 
        lwd = 2, xlab = expression(theta), ylab = expression(phi))
points(to_rad(samp), pch = 16)