Nothing
tests/testthat/tests_bwd.R
In polykde: Polyspherical Kernel Density Estimation
## Constants kernel
d <- sample(1:4, size = 1, replace = TRUE)
test_that("b_d(L) equals definition for kernel = 1, 2", {
expect_equal(b_d(d = d, kernel = 1),
DirStats::b_L(L = function(r) exp(-r), q = d) / d,
tolerance = 1e-6)
expect_equal(b_d(d = d, kernel = 2),
DirStats::b_L(L = function(r) (1 - r) * (r <= 1), q = d) / d,
tolerance = 1e-6)
})
test_that("v_d(L) equals definition for kernel = 1, 2", {
expect_equal(v_d(d = d, kernel = 1),
DirStats::d_L(L = function(r) exp(-r), q = d) /
DirStats::lambda_L(L = function(r) exp(-r), q = d),
tolerance = 1e-6)
expect_equal(v_d(d = d, kernel = 2),
DirStats::d_L(L = function(r) (1 - r) * (r <= 1), q = d) /
DirStats::lambda_L(L = function(r) (1 - r) * (r <= 1), q = d),
tolerance = 1e-6)
})
## CV bandwidth selectors
seed <- 30
set.seed(seed, kind = "Mersenne-Twister")
r <- 2
d <- sample(1:3, size = r, replace = TRUE)
h <- sample(c(0.25, 0.5, 0.75), size = r, replace = TRUE)
n <- 20
mu <- r_unif_polysph(n = 5, d = d)
X <- r_kde_polysph(n = n, X = mu, d = d, h = h)
# CV helper functions
M <- 1e4
seed <- 30
set.seed(seed, kind = "Mersenne-Twister")
mc_samp <- r_unif_polysph(n = M, d = d)
cv_naive <- function(h, X, d, mc_samp, kde_samp = FALSE) {
if (kde_samp) {
set.seed(seed, kind = "Mersenne-Twister")
mc_kde_samp <- r_kde_polysph(n = M, X = X, d = d, h = h, kernel = 1)
cv_1 <- mean(kde_polysph(x = mc_kde_samp, X = X, d = d, kernel = 1, h = h,
wrt_unif = FALSE))
cv_2 <- 2 * mean(exp(log_cv_kde_polysph(X = X, d = d, kernel = 1, h = h,
wrt_unif = FALSE)))
} else {
cv_1 <- prod(rotasym::w_p(p = d + 1)) *
mean(kde_polysph(x = mc_samp, X = X, d = d, kernel = 1, h = h,
wrt_unif = FALSE)^2)
cv_2 <- 2 * mean(exp(log_cv_kde_polysph(X = X, d = d, kernel = 1, h = h,
wrt_unif = FALSE)))
}
return(cv_1 - cv_2)
}
# cv_naive_curve <- function(h1, kde_samp = FALSE) sapply(h1, function(hh1) {
# cv_naive(h = rep(hh1, r), X = X, d = d, mc_samp = mc_samp,
# kde_samp = kde_samp)
# })
# cv_bw_cv_curve <- function(h1) sapply(h1, function(hh1) {
# bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
# bw0 = rep(hh1, r), M = M, control = list(maxit = 0),
# method = "BFGS", exact_vmf = FALSE, seed_mc = seed)$opt$value
# })
# cv_bw_cv_vmf_curve <- function(h1) sapply(h1, function(hh1) {
# bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
# bw0 = rep(hh1, r), control = list(maxit = 0),
# method = "BFGS", exact_vmf = TRUE)$opt$value
# })
#
# # Visualization of LSCV functions
# curve(cv_naive_curve(x, kde_samp = FALSE), from = 0.2, to = 1, n = 100,
# ylab = "CV loss")
# curve(cv_naive_curve(x, kde_samp = TRUE), from = 0.2, to = 1, n = 100,
# add = TRUE, lty = 2)
# curve(cv_bw_cv_curve, from = 0.2, to = 1, n = 100, add = TRUE, col = 2)
# curve(cv_bw_cv_vmf_curve, from = 0.2, to = 1, n = 100, add = TRUE, col = 3)
test_that("bw_cv_polysph(type = \"LCV\") in sequential and parallel mode", {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV",
control = list(maxit = 1e3))$opt$value,
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV",
control = list(maxit = 1e3), ncores = 2)$opt$value,
tolerance = 1e-2)
})
test_that("bw_cv_polysph(type = \"LSCV\") in sequential and parallel mode", {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
control = list(maxit = 1e3), seed_mc = 1)$opt$value,
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
control = list(maxit = 1e3), seed_mc = 1,
ncores = 2)$opt$value,
tolerance = 1e-2)
})
test_that("bw_cv_polysph(type = \"LSCV\", imp_mc = TRUE) loss", {
for (f in c(0.25, 0.5, 1, 2)) {
expect_equal(
cv_naive(h = f * h, X = X, d = d, kde_samp = TRUE),
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
bw0 = f * h, M = M, control = list(maxit = 0),
method = "BFGS", exact_vmf = FALSE, imp_mc = TRUE,
seed_mc = seed)$opt$value)
}
})
test_that("bw_cv_polysph(type = \"LSCV\", imp_mc = FALSE) loss", {
for (f in c(0.25, 0.5, 1, 2)) {
expect_equal(
cv_naive(h = f * h, X = X, d = d, mc_samp = mc_samp, kde_samp = FALSE),
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
bw0 = f * h, M = M, control = list(maxit = 0),
method = "BFGS", exact_vmf = FALSE, imp_mc = FALSE,
seed_mc = seed)$opt$value,
tolerance = 5e-2)
}
})
test_that("bw_cv_polysph(type = \"LSCV\", exact_vmf = TRUE) loss", {
for (f in c(0.25, 0.5, 1, 2)) {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
bw0 = f * h, M = M, control = list(maxit = 0),
method = "BFGS", exact_vmf = FALSE,
seed_mc = seed)$opt$value,
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
bw0 = f * h, M = M, control = list(maxit = 0),
method = "BFGS", exact_vmf = TRUE)$opt$value,
tolerance = 5e-2)
}
})
test_that("bw_cv_polysph(type = \"LCV\", common_h = TRUE)", {
for (f in c(0.25, 0.5, 1, 2)) {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV",
bw0 = f * h, M = M, control = list(maxit = 0),
method = "BFGS", common_h = TRUE)$opt$value,
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV",
bw0 = f * rep(mean(h), r), M = M, control = list(maxit = 0),
method = "BFGS", common_h = FALSE)$opt$value,
tolerance = 5e-2)
}
})
test_that("bw_cv_polysph(type = \"LSCV\", common_h = TRUE)", {
for (f in c(0.25, 0.5, 1, 2)) {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
bw0 = f * h, M = M, control = list(maxit = 0),
method = "BFGS", exact_vmf = TRUE, common_h = TRUE,
seed_mc = seed)$opt$value,
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
bw0 = f * rep(mean(h), r), M = M, control = list(maxit = 0),
method = "BFGS", exact_vmf = TRUE,
common_h = FALSE)$opt$value,
tolerance = 5e-2)
}
})
test_that("bw_cv_polysph() with bw0 vector and bw0 matrix", {
bw0_vec_1 <- h
bw0_vec_2 <- 0.1 * h
bw0_vec_3 <- 5 * h
bw0_mat <- rbind(bw0_vec_1, bw0_vec_2, bw0_vec_3)
expect_true(
bw_cv_polysph(X = X, d = d, bw0 = bw0_vec_1,
control = list(maxit = 0))$opt$value >=
bw_cv_polysph(X = X, d = d, bw0 = bw0_mat,
control = list(maxit = 0))$opt$value)
expect_true(
bw_cv_polysph(X = X, d = d, bw0 = bw0_vec_2,
control = list(maxit = 0))$opt$value >=
bw_cv_polysph(X = X, d = d, bw0 = bw0_mat,
control = list(maxit = 0))$opt$value)
expect_true(
bw_cv_polysph(X = X, d = d, bw0 = bw0_vec_3,
control = list(maxit = 0))$opt$value >=
bw_cv_polysph(X = X, d = d, bw0 = bw0_mat,
control = list(maxit = 0))$opt$value)
})
test_that("bw_cv_polysph() with bw0 vector and bw0 matrix with common_h", {
bw0_vec_1 <- h
bw0_vec_2 <- 0.1 * h
bw0_vec_3 <- 5 * h
bw0_mat <- rbind(bw0_vec_1, bw0_vec_2, bw0_vec_3)
suppressWarnings({
expect_true(
bw_cv_polysph(X = X, d = d, bw0 = bw0_vec_1, control = list(maxit = 0),
common_h = TRUE)$opt$value >=
bw_cv_polysph(X = X, d = d, bw0 = bw0_mat, control = list(maxit = 0),
common_h = TRUE)$opt$value)
expect_true(
bw_cv_polysph(X = X, d = d, bw0 = bw0_vec_2, control = list(maxit = 0),
common_h = TRUE)$opt$value >=
bw_cv_polysph(X = X, d = d, bw0 = bw0_mat, control = list(maxit = 0),
common_h = TRUE)$opt$value)
expect_true(
bw_cv_polysph(X = X, d = d, bw0 = bw0_vec_3, control = list(maxit = 0),
common_h = TRUE)$opt$value >=
bw_cv_polysph(X = X, d = d, bw0 = bw0_mat, control = list(maxit = 0),
common_h = TRUE)$opt$value)
})
})
test_that("bw_cv_polysph(type = \"LCV\") with optim/nlm", {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV", opt = "nlm")$bw,
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV", opt = "optim")$bw,
tolerance = 1e-2)
expect_error(bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV",
opt = "nlm", ncores = 2))
})
## CV equivalence with DirStats::bw_dir_*cv()
r <- 1
d <- 2
h <- 0.25
n <- 20
mu <- r_unif_polysph(n = 5, d = d)
X <- r_kde_polysph(n = n, X = mu, d = d, h = h)
test_that("bw_cv_polysph(type = \"LCV\") equals DirStats::bw_dir_lcv()", {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV",
method = "L-BFGS-B", bw0 = h)$bw,
DirStats::bw_dir_lcv(data = X, optim = TRUE, optim_par = h)$h_opt,
tolerance = 1e-3)
})
test_that("bw_cv_polysph(type = \"LSCV\") equals DirStats::bw_dir_lscv()", {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
method = "L-BFGS-B", bw0 = 0.5, exact_vmf = TRUE)$bw,
DirStats::bw_dir_lscv(data = X, optim = TRUE, optim_par = 0.5)$h_opt,
tolerance = 5e-2)
})
## Curvature matrices for von Mises--Fisher
r <- rpois(1, lambda = 3) + 1
d <- rpois(r, lambda = 2) + 1
mu <- r_unif_polysph(n = 1, d = d)
kappa <- r * runif(r)
ind_dj <- comp_ind_dj(d = d)
# Approximate curvature using weighted inverse sampling
# \int_{S^{d_1, ..., d_r}} t(x)t(x)' dx
# = \int \int_{S^{d_1, ..., d_r}} [t(x) / f(x)][t(x) / f(x)]' f(x)^2 dx
# ≈ 1 / M \sum_{i = 1}^M [t(X_i) / f(X_i)][t(X_i) / f(X_i)]' * f(X_i)
# where t(x) = (tr(H_11 f(x)), ..., tr(H_rr f(x)))'.
# t() for vMF
hess_vmf <- function(x) grad_hess_kde_polysph(x = rbind(x), X = mu, d = d,
h = 1 / sqrt(kappa),
norm_grad_hess = TRUE)$hess[1, , ]
t_vmf <- function(x) {
H <- hess_vmf(x)
sapply(1:r, function(j) {
ind_j <- (ind_dj[j] + 1):ind_dj[j + 1]
sum(diag(H[ind_j, ind_j]))
})
}
# Integral estimation
N <- 1e4
X_vmf <- r_vmf_polysph(n = N, d = d, mu = mu, kappa = kappa)
f_X_vmf <- drop(kde_polysph(x = X_vmf, X = mu, d = d, h = 1 / sqrt(kappa)))
tt_X_vmf <- t(rbind(apply(X_vmf, 1, function(xi) tcrossprod(t_vmf(x = xi)))))
R_X_vmf <- matrix(colMeans(tt_X_vmf * f_X_vmf), nrow = r, ncol = r)
test_that("curv_vmf_polysph() curvature matrix against its numerical version", {
expect_equal(R_X_vmf,
curv_vmf_polysph(kappa = kappa, d = d),
tolerance = 1e-2)
})
test_that("sum(curv_vmf_polysph()) curvature term against its numerical
version", {
expect_equal(sum(R_X_vmf),
sum(curv_vmf_polysph(kappa = kappa, d = d)),
tolerance = 5e-2)
})
test_that("curv_vmf_polysph() is coherent with scalar/vector d", {
expect_equal(curv_vmf_polysph(kappa = kappa, d = rep(d[1], r)),
curv_vmf_polysph(kappa = kappa, d = d[1]))
})
test_that("curv_vmf_polysph() equivalence with DirStats::R_Psi_mixvmf()", {
expect_equal(drop(curv_vmf_polysph(kappa = kappa[1], d = d[1])),
d[1]^2 * DirStats::R_Psi_mixvmf(q = d[1],
mu = rbind(c(rep(0, d[1]), 1)),
kappa = kappa[1], p = 1))
})
## Objective functions for plug-in bandwidth selectors
# Setup
r <- rpois(1, lambda = 5) + 1
d <- rpois(r, lambda = 2) + 1
h <- runif(r)
n <- 50
kappa <- 1:r
mu <- r_unif_polysph(n = 5, d = d)
X <- r_kde_polysph(n = n, X = mu, d = d, h = h)
ind <- cumsum(c(1, d + 1))
# Common objects
R_kappa <- curv_vmf_polysph(kappa = kappa, d = d)
bias2 <- tcrossprod(b_d(kernel = 1, d = d)) * R_kappa
var <- prod(v_d(kernel = 1, d = d)) / n
# AMISE and derivative
log_amise <- function(h) {
var2 <- var / abs(prod(h^d))
obj <- sum(tcrossprod(h^2) * bias2) + var2
attr(obj, "gradient") <- (bias2 %*% (4 * h^2) * h - d * var2 / abs(h)) / obj
return(log(obj))
}
# Common objects
log_R_kappa <- curv_vmf_polysph(kappa = kappa, d = d, log = TRUE)
log_bias2 <- log(tcrossprod(b_d(kernel = 1, d = d))) + log_R_kappa
log_var <- sum(log(v_d(kernel = 1, d = d))) - log(n)
# log(exp(log_x) + exp(y))
log_sum_exp <- function(x) {
M <- max(x)
M + log(sum(exp(x - M)))
}
# AMISE and derivative
log_amise_stable <- function(h) {
log_h <- log(abs(h))
log_var2 <- log_var - sum(d * log_h)
logs <- log(tcrossprod(h^2)) + log_bias2 - log_var2
log_obj <- log_var2 + polykde:::log_sum_exp(logs = c(logs, 0))
attr(log_obj, "gradient") <-
exp(log(4) + log_bias2 - log_obj + log_h) %*% h^2 -
exp(log(d) + log_var2 - log_obj - log_h)
return(log_obj)
}
log_amise_stable_log_h <- function(log_h) {
h <- exp(log_h)
log_var2 <- log_var - sum(d * log_h)
logs <- log(tcrossprod(h^2)) + log_bias2 - log_var2
log_obj <- log_var2 + polykde:::log_sum_exp(logs = c(logs, 0))
attr(log_obj, "gradient") <-
(exp(log(4) + log_bias2 - log_obj + log_h) %*% h^2 -
exp(log(d) + log_var2 - log_obj - log_h)) * h
return(log_obj)
}
test_that("Coherence between bw_rot_polysph() and bw_mrot_polysph()", {
expect_equal(
sapply(seq_along(d), function(j) {
data <- X[, ind[j]:(ind[j + 1] - 1)]
bw_rot_polysph(X = data, d = d[j])$bw
}),
bw_mrot_polysph(X = X, d = d),
tolerance = 1e-6)
})
test_that("Coherence between bw_rot_polysph() and DirStats::bw_dir_rot()", {
expect_equal(
sapply(seq_along(d), function(j) {
data <- X[, ind[j]:(ind[j + 1] - 1)]
bw_rot_polysph(X = data, d = d[j])$bw
}),
sapply(seq_along(d), function(j) {
data <- X[, ind[j]:(ind[j + 1] - 1)]
DirStats::bw_dir_rot(data = data)
}), tolerance = 1e-4)
})
test_that("Derivatives of log_amise_*", {
for (i in 1:5) {
h <- 2 * runif(r)
expect_equal(
drop(attr(log_amise(h), which = "gradient")),
numDeriv::grad(func = log_amise, x = h, method = "Richardson",
method.args = list(eps = 1e-10)),
tolerance = 5e-6
)
expect_equal(
drop(attr(log_amise_stable(h), which = "gradient")),
numDeriv::grad(func = log_amise_stable, x = h, method = "Richardson",
method.args = list(eps = 1e-10)),
tolerance = 5e-6
)
expect_equal(
drop(attr(log_amise_stable_log_h(log(h)), which = "gradient")),
numDeriv::grad(func = log_amise_stable_log_h, x = log(h),
method = "Richardson", method.args = list(eps = 1e-10)),
tolerance = 5e-6
)
}
})
test_that("log_amise vs. log_amise_stable vs. log_amise_stable_log_h", {
for (i in 1:5) {
h <- 2 * runif(r)
expect_equal(
log_amise(h),
log_amise_stable(h)
)
expect_equal(
log_amise(h),
{lh <- log_amise_stable_log_h(log(h))
attr(lh, "gradient") <- attr(lh, "gradient") / h
lh}
)
}
})
test_that("log_amise vs. log_amise_stable", {
for (i in 1:5) {
h <- 2 * runif(r)
expect_equal(
log_amise(h),
log_amise_stable(h)
)
}
})
test_that("Same optimization with log_amise_stable or log_amise_stable_log_h", {
h0 <- 2 * runif(r)
opt1 <- nlm(f = log_amise_stable_log_h, p = log(h0))
opt2 <- nlm(f = log_amise_stable, p = exp(opt1$estimate))
opt3 <- nlm(f = log_amise_stable, p = h0)
expect_true(
(abs(opt1$minimum - opt2$minimum) < 1e-9) ||
(max(abs(exp(opt1$est) - abs(opt2$est))) < 1e-9)
)
expect_lt(opt1$minimum - 1e-9, opt3$minimum)
})
test_that("bw_rot_polysph() with bw0 vector and bw0 matrix", {
bw0_vec_1 <- h
bw0_vec_2 <- 0.1 * h
bw0_vec_3 <- 5 * h
bw0_mat <- rbind(bw0_vec_1, bw0_vec_2, bw0_vec_3)
expect_true(bw_rot_polysph(X = X, d = d, bw0 = bw0_vec_1,
iterlim = 1)$opt$minimum + 0.1 >=
bw_rot_polysph(X = X, d = d, bw0 = bw0_mat,
iterlim = 1)$opt$minimum)
expect_true(bw_rot_polysph(X = X, d = d, bw0 = bw0_vec_2,
iterlim = 1)$opt$minimum + 0.1 >=
bw_rot_polysph(X = X, d = d, bw0 = bw0_mat,
iterlim = 1)$opt$minimum)
expect_true(bw_rot_polysph(X = X, d = d, bw0 = bw0_vec_3,
iterlim = 1)$opt$minimum + 0.1 >=
bw_rot_polysph(X = X, d = d, bw0 = bw0_mat,
iterlim = 1)$opt$minimum)
})
## Plug-in bandwidth selectors
test_that("Same result for vMF kernel with kernel_type = 1,2", {
# Setup
r <- rpois(1, lambda = 5) + 1
d <- rpois(r, lambda = 2) + 1
h <- runif(r)
n <- 50
kappa <- 1:r
mu <- r_unif_polysph(n = 5, d = d)
X <- r_kde_polysph(n = n, X = mu, d = d, h = h)
expect_equal(
bw_rot_polysph(X = X, d = d, kernel = 1, kernel_type = 1)$bw,
bw_rot_polysph(X = X, d = d, kernel = 1, kernel_type = 2)$bw
)
})
test_that("Same result for Epa and sfp kernel with kernel_type = 1,2
and r = 1", {
# Setup
d <- rpois(1, lambda = 2) + 1
mu <- r_unif_polysph(n = 1, d = d)
X <- r_kde_polysph(n = 50, X = mu, d = d, h = runif(1))
expect_equal(
bw_rot_polysph(X = X, d = d, kernel = 2, kernel_type = 1)$bw,
bw_rot_polysph(X = X, d = d, kernel = 2, kernel_type = 2)$bw
)
expect_equal(
bw_rot_polysph(X = X, d = d, kernel = 3, kernel_type = 1)$bw,
bw_rot_polysph(X = X, d = d, kernel = 3, kernel_type = 2)$bw
)
})
test_that("Same result with kappa precomputed", {
# Setup
r <- 2
d <- rpois(r, lambda = 2) + 1
mu <- r_unif_polysph(n = 1, d = d)
X <- r_kde_polysph(n = 50, X = mu, d = d, h = runif(r))
ind <- cumsum(c(1, d + 1))
kappa <- sapply(seq_along(d), function(j) {
# Prepare data + fit vMF
data <- X[, ind[j]:(ind[j + 1] - 1)]
min(DirStats::norm2(movMF::movMF(x = data, k = 1, type = "S",
maxit = 300)$theta),
5e4) # Breaking point for later Bessels
})
expect_equal(
bw_rot_polysph(X = X + 1, d = d, kappa = kappa)$bw,
bw_rot_polysph(X = X, d = d, kappa = NULL)$bw
)
expect_equal(
bw_mrot_polysph(X = X + 1, d = d, kappa = kappa),
bw_mrot_polysph(X = X, d = d, kappa = NULL)
)
})
Try the polykde package in your browser
Any scripts or data that you put into this service are public.
polykde documentation built on April 16, 2025, 1:11 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## Constants kernel
d <- sample(1:4, size = 1, replace = TRUE)
test_that("b_d(L) equals definition for kernel = 1, 2", {
expect_equal(b_d(d = d, kernel = 1),
DirStats::b_L(L = function(r) exp(-r), q = d) / d,
tolerance = 1e-6)
expect_equal(b_d(d = d, kernel = 2),
DirStats::b_L(L = function(r) (1 - r) * (r <= 1), q = d) / d,
tolerance = 1e-6)
})
test_that("v_d(L) equals definition for kernel = 1, 2", {
expect_equal(v_d(d = d, kernel = 1),
DirStats::d_L(L = function(r) exp(-r), q = d) /
DirStats::lambda_L(L = function(r) exp(-r), q = d),
tolerance = 1e-6)
expect_equal(v_d(d = d, kernel = 2),
DirStats::d_L(L = function(r) (1 - r) * (r <= 1), q = d) /
DirStats::lambda_L(L = function(r) (1 - r) * (r <= 1), q = d),
tolerance = 1e-6)
})
## CV bandwidth selectors
seed <- 30
set.seed(seed, kind = "Mersenne-Twister")
r <- 2
d <- sample(1:3, size = r, replace = TRUE)
h <- sample(c(0.25, 0.5, 0.75), size = r, replace = TRUE)
n <- 20
mu <- r_unif_polysph(n = 5, d = d)
X <- r_kde_polysph(n = n, X = mu, d = d, h = h)
# CV helper functions
M <- 1e4
seed <- 30
set.seed(seed, kind = "Mersenne-Twister")
mc_samp <- r_unif_polysph(n = M, d = d)
cv_naive <- function(h, X, d, mc_samp, kde_samp = FALSE) {
if (kde_samp) {
set.seed(seed, kind = "Mersenne-Twister")
mc_kde_samp <- r_kde_polysph(n = M, X = X, d = d, h = h, kernel = 1)
cv_1 <- mean(kde_polysph(x = mc_kde_samp, X = X, d = d, kernel = 1, h = h,
wrt_unif = FALSE))
cv_2 <- 2 * mean(exp(log_cv_kde_polysph(X = X, d = d, kernel = 1, h = h,
wrt_unif = FALSE)))
} else {
cv_1 <- prod(rotasym::w_p(p = d + 1)) *
mean(kde_polysph(x = mc_samp, X = X, d = d, kernel = 1, h = h,
wrt_unif = FALSE)^2)
cv_2 <- 2 * mean(exp(log_cv_kde_polysph(X = X, d = d, kernel = 1, h = h,
wrt_unif = FALSE)))
}
return(cv_1 - cv_2)
}
# cv_naive_curve <- function(h1, kde_samp = FALSE) sapply(h1, function(hh1) {
# cv_naive(h = rep(hh1, r), X = X, d = d, mc_samp = mc_samp,
# kde_samp = kde_samp)
# })
# cv_bw_cv_curve <- function(h1) sapply(h1, function(hh1) {
# bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
# bw0 = rep(hh1, r), M = M, control = list(maxit = 0),
# method = "BFGS", exact_vmf = FALSE, seed_mc = seed)$opt$value
# })
# cv_bw_cv_vmf_curve <- function(h1) sapply(h1, function(hh1) {
# bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
# bw0 = rep(hh1, r), control = list(maxit = 0),
# method = "BFGS", exact_vmf = TRUE)$opt$value
# })
#
# # Visualization of LSCV functions
# curve(cv_naive_curve(x, kde_samp = FALSE), from = 0.2, to = 1, n = 100,
# ylab = "CV loss")
# curve(cv_naive_curve(x, kde_samp = TRUE), from = 0.2, to = 1, n = 100,
# add = TRUE, lty = 2)
# curve(cv_bw_cv_curve, from = 0.2, to = 1, n = 100, add = TRUE, col = 2)
# curve(cv_bw_cv_vmf_curve, from = 0.2, to = 1, n = 100, add = TRUE, col = 3)
test_that("bw_cv_polysph(type = \"LCV\") in sequential and parallel mode", {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV",
control = list(maxit = 1e3))$opt$value,
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV",
control = list(maxit = 1e3), ncores = 2)$opt$value,
tolerance = 1e-2)
})
test_that("bw_cv_polysph(type = \"LSCV\") in sequential and parallel mode", {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
control = list(maxit = 1e3), seed_mc = 1)$opt$value,
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
control = list(maxit = 1e3), seed_mc = 1,
ncores = 2)$opt$value,
tolerance = 1e-2)
})
test_that("bw_cv_polysph(type = \"LSCV\", imp_mc = TRUE) loss", {
for (f in c(0.25, 0.5, 1, 2)) {
expect_equal(
cv_naive(h = f * h, X = X, d = d, kde_samp = TRUE),
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
bw0 = f * h, M = M, control = list(maxit = 0),
method = "BFGS", exact_vmf = FALSE, imp_mc = TRUE,
seed_mc = seed)$opt$value)
}
})
test_that("bw_cv_polysph(type = \"LSCV\", imp_mc = FALSE) loss", {
for (f in c(0.25, 0.5, 1, 2)) {
expect_equal(
cv_naive(h = f * h, X = X, d = d, mc_samp = mc_samp, kde_samp = FALSE),
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
bw0 = f * h, M = M, control = list(maxit = 0),
method = "BFGS", exact_vmf = FALSE, imp_mc = FALSE,
seed_mc = seed)$opt$value,
tolerance = 5e-2)
}
})
test_that("bw_cv_polysph(type = \"LSCV\", exact_vmf = TRUE) loss", {
for (f in c(0.25, 0.5, 1, 2)) {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
bw0 = f * h, M = M, control = list(maxit = 0),
method = "BFGS", exact_vmf = FALSE,
seed_mc = seed)$opt$value,
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
bw0 = f * h, M = M, control = list(maxit = 0),
method = "BFGS", exact_vmf = TRUE)$opt$value,
tolerance = 5e-2)
}
})
test_that("bw_cv_polysph(type = \"LCV\", common_h = TRUE)", {
for (f in c(0.25, 0.5, 1, 2)) {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV",
bw0 = f * h, M = M, control = list(maxit = 0),
method = "BFGS", common_h = TRUE)$opt$value,
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV",
bw0 = f * rep(mean(h), r), M = M, control = list(maxit = 0),
method = "BFGS", common_h = FALSE)$opt$value,
tolerance = 5e-2)
}
})
test_that("bw_cv_polysph(type = \"LSCV\", common_h = TRUE)", {
for (f in c(0.25, 0.5, 1, 2)) {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
bw0 = f * h, M = M, control = list(maxit = 0),
method = "BFGS", exact_vmf = TRUE, common_h = TRUE,
seed_mc = seed)$opt$value,
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
bw0 = f * rep(mean(h), r), M = M, control = list(maxit = 0),
method = "BFGS", exact_vmf = TRUE,
common_h = FALSE)$opt$value,
tolerance = 5e-2)
}
})
test_that("bw_cv_polysph() with bw0 vector and bw0 matrix", {
bw0_vec_1 <- h
bw0_vec_2 <- 0.1 * h
bw0_vec_3 <- 5 * h
bw0_mat <- rbind(bw0_vec_1, bw0_vec_2, bw0_vec_3)
expect_true(
bw_cv_polysph(X = X, d = d, bw0 = bw0_vec_1,
control = list(maxit = 0))$opt$value >=
bw_cv_polysph(X = X, d = d, bw0 = bw0_mat,
control = list(maxit = 0))$opt$value)
expect_true(
bw_cv_polysph(X = X, d = d, bw0 = bw0_vec_2,
control = list(maxit = 0))$opt$value >=
bw_cv_polysph(X = X, d = d, bw0 = bw0_mat,
control = list(maxit = 0))$opt$value)
expect_true(
bw_cv_polysph(X = X, d = d, bw0 = bw0_vec_3,
control = list(maxit = 0))$opt$value >=
bw_cv_polysph(X = X, d = d, bw0 = bw0_mat,
control = list(maxit = 0))$opt$value)
})
test_that("bw_cv_polysph() with bw0 vector and bw0 matrix with common_h", {
bw0_vec_1 <- h
bw0_vec_2 <- 0.1 * h
bw0_vec_3 <- 5 * h
bw0_mat <- rbind(bw0_vec_1, bw0_vec_2, bw0_vec_3)
suppressWarnings({
expect_true(
bw_cv_polysph(X = X, d = d, bw0 = bw0_vec_1, control = list(maxit = 0),
common_h = TRUE)$opt$value >=
bw_cv_polysph(X = X, d = d, bw0 = bw0_mat, control = list(maxit = 0),
common_h = TRUE)$opt$value)
expect_true(
bw_cv_polysph(X = X, d = d, bw0 = bw0_vec_2, control = list(maxit = 0),
common_h = TRUE)$opt$value >=
bw_cv_polysph(X = X, d = d, bw0 = bw0_mat, control = list(maxit = 0),
common_h = TRUE)$opt$value)
expect_true(
bw_cv_polysph(X = X, d = d, bw0 = bw0_vec_3, control = list(maxit = 0),
common_h = TRUE)$opt$value >=
bw_cv_polysph(X = X, d = d, bw0 = bw0_mat, control = list(maxit = 0),
common_h = TRUE)$opt$value)
})
})
test_that("bw_cv_polysph(type = \"LCV\") with optim/nlm", {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV", opt = "nlm")$bw,
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV", opt = "optim")$bw,
tolerance = 1e-2)
expect_error(bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV",
opt = "nlm", ncores = 2))
})
## CV equivalence with DirStats::bw_dir_*cv()
r <- 1
d <- 2
h <- 0.25
n <- 20
mu <- r_unif_polysph(n = 5, d = d)
X <- r_kde_polysph(n = n, X = mu, d = d, h = h)
test_that("bw_cv_polysph(type = \"LCV\") equals DirStats::bw_dir_lcv()", {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LCV",
method = "L-BFGS-B", bw0 = h)$bw,
DirStats::bw_dir_lcv(data = X, optim = TRUE, optim_par = h)$h_opt,
tolerance = 1e-3)
})
test_that("bw_cv_polysph(type = \"LSCV\") equals DirStats::bw_dir_lscv()", {
expect_equal(
bw_cv_polysph(X = X, d = d, kernel = 1, type = "LSCV",
method = "L-BFGS-B", bw0 = 0.5, exact_vmf = TRUE)$bw,
DirStats::bw_dir_lscv(data = X, optim = TRUE, optim_par = 0.5)$h_opt,
tolerance = 5e-2)
})
## Curvature matrices for von Mises--Fisher
r <- rpois(1, lambda = 3) + 1
d <- rpois(r, lambda = 2) + 1
mu <- r_unif_polysph(n = 1, d = d)
kappa <- r * runif(r)
ind_dj <- comp_ind_dj(d = d)
# Approximate curvature using weighted inverse sampling
# \int_{S^{d_1, ..., d_r}} t(x)t(x)' dx
# = \int \int_{S^{d_1, ..., d_r}} [t(x) / f(x)][t(x) / f(x)]' f(x)^2 dx
# ≈ 1 / M \sum_{i = 1}^M [t(X_i) / f(X_i)][t(X_i) / f(X_i)]' * f(X_i)
# where t(x) = (tr(H_11 f(x)), ..., tr(H_rr f(x)))'.
# t() for vMF
hess_vmf <- function(x) grad_hess_kde_polysph(x = rbind(x), X = mu, d = d,
h = 1 / sqrt(kappa),
norm_grad_hess = TRUE)$hess[1, , ]
t_vmf <- function(x) {
H <- hess_vmf(x)
sapply(1:r, function(j) {
ind_j <- (ind_dj[j] + 1):ind_dj[j + 1]
sum(diag(H[ind_j, ind_j]))
})
}
# Integral estimation
N <- 1e4
X_vmf <- r_vmf_polysph(n = N, d = d, mu = mu, kappa = kappa)
f_X_vmf <- drop(kde_polysph(x = X_vmf, X = mu, d = d, h = 1 / sqrt(kappa)))
tt_X_vmf <- t(rbind(apply(X_vmf, 1, function(xi) tcrossprod(t_vmf(x = xi)))))
R_X_vmf <- matrix(colMeans(tt_X_vmf * f_X_vmf), nrow = r, ncol = r)
test_that("curv_vmf_polysph() curvature matrix against its numerical version", {
expect_equal(R_X_vmf,
curv_vmf_polysph(kappa = kappa, d = d),
tolerance = 1e-2)
})
test_that("sum(curv_vmf_polysph()) curvature term against its numerical
version", {
expect_equal(sum(R_X_vmf),
sum(curv_vmf_polysph(kappa = kappa, d = d)),
tolerance = 5e-2)
})
test_that("curv_vmf_polysph() is coherent with scalar/vector d", {
expect_equal(curv_vmf_polysph(kappa = kappa, d = rep(d[1], r)),
curv_vmf_polysph(kappa = kappa, d = d[1]))
})
test_that("curv_vmf_polysph() equivalence with DirStats::R_Psi_mixvmf()", {
expect_equal(drop(curv_vmf_polysph(kappa = kappa[1], d = d[1])),
d[1]^2 * DirStats::R_Psi_mixvmf(q = d[1],
mu = rbind(c(rep(0, d[1]), 1)),
kappa = kappa[1], p = 1))
})
## Objective functions for plug-in bandwidth selectors
# Setup
r <- rpois(1, lambda = 5) + 1
d <- rpois(r, lambda = 2) + 1
h <- runif(r)
n <- 50
kappa <- 1:r
mu <- r_unif_polysph(n = 5, d = d)
X <- r_kde_polysph(n = n, X = mu, d = d, h = h)
ind <- cumsum(c(1, d + 1))
# Common objects
R_kappa <- curv_vmf_polysph(kappa = kappa, d = d)
bias2 <- tcrossprod(b_d(kernel = 1, d = d)) * R_kappa
var <- prod(v_d(kernel = 1, d = d)) / n
# AMISE and derivative
log_amise <- function(h) {
var2 <- var / abs(prod(h^d))
obj <- sum(tcrossprod(h^2) * bias2) + var2
attr(obj, "gradient") <- (bias2 %*% (4 * h^2) * h - d * var2 / abs(h)) / obj
return(log(obj))
}
# Common objects
log_R_kappa <- curv_vmf_polysph(kappa = kappa, d = d, log = TRUE)
log_bias2 <- log(tcrossprod(b_d(kernel = 1, d = d))) + log_R_kappa
log_var <- sum(log(v_d(kernel = 1, d = d))) - log(n)
# log(exp(log_x) + exp(y))
log_sum_exp <- function(x) {
M <- max(x)
M + log(sum(exp(x - M)))
}
# AMISE and derivative
log_amise_stable <- function(h) {
log_h <- log(abs(h))
log_var2 <- log_var - sum(d * log_h)
logs <- log(tcrossprod(h^2)) + log_bias2 - log_var2
log_obj <- log_var2 + polykde:::log_sum_exp(logs = c(logs, 0))
attr(log_obj, "gradient") <-
exp(log(4) + log_bias2 - log_obj + log_h) %*% h^2 -
exp(log(d) + log_var2 - log_obj - log_h)
return(log_obj)
}
log_amise_stable_log_h <- function(log_h) {
h <- exp(log_h)
log_var2 <- log_var - sum(d * log_h)
logs <- log(tcrossprod(h^2)) + log_bias2 - log_var2
log_obj <- log_var2 + polykde:::log_sum_exp(logs = c(logs, 0))
attr(log_obj, "gradient") <-
(exp(log(4) + log_bias2 - log_obj + log_h) %*% h^2 -
exp(log(d) + log_var2 - log_obj - log_h)) * h
return(log_obj)
}
test_that("Coherence between bw_rot_polysph() and bw_mrot_polysph()", {
expect_equal(
sapply(seq_along(d), function(j) {
data <- X[, ind[j]:(ind[j + 1] - 1)]
bw_rot_polysph(X = data, d = d[j])$bw
}),
bw_mrot_polysph(X = X, d = d),
tolerance = 1e-6)
})
test_that("Coherence between bw_rot_polysph() and DirStats::bw_dir_rot()", {
expect_equal(
sapply(seq_along(d), function(j) {
data <- X[, ind[j]:(ind[j + 1] - 1)]
bw_rot_polysph(X = data, d = d[j])$bw
}),
sapply(seq_along(d), function(j) {
data <- X[, ind[j]:(ind[j + 1] - 1)]
DirStats::bw_dir_rot(data = data)
}), tolerance = 1e-4)
})
test_that("Derivatives of log_amise_*", {
for (i in 1:5) {
h <- 2 * runif(r)
expect_equal(
drop(attr(log_amise(h), which = "gradient")),
numDeriv::grad(func = log_amise, x = h, method = "Richardson",
method.args = list(eps = 1e-10)),
tolerance = 5e-6
)
expect_equal(
drop(attr(log_amise_stable(h), which = "gradient")),
numDeriv::grad(func = log_amise_stable, x = h, method = "Richardson",
method.args = list(eps = 1e-10)),
tolerance = 5e-6
)
expect_equal(
drop(attr(log_amise_stable_log_h(log(h)), which = "gradient")),
numDeriv::grad(func = log_amise_stable_log_h, x = log(h),
method = "Richardson", method.args = list(eps = 1e-10)),
tolerance = 5e-6
)
}
})
test_that("log_amise vs. log_amise_stable vs. log_amise_stable_log_h", {
for (i in 1:5) {
h <- 2 * runif(r)
expect_equal(
log_amise(h),
log_amise_stable(h)
)
expect_equal(
log_amise(h),
{lh <- log_amise_stable_log_h(log(h))
attr(lh, "gradient") <- attr(lh, "gradient") / h
lh}
)
}
})
test_that("log_amise vs. log_amise_stable", {
for (i in 1:5) {
h <- 2 * runif(r)
expect_equal(
log_amise(h),
log_amise_stable(h)
)
}
})
test_that("Same optimization with log_amise_stable or log_amise_stable_log_h", {
h0 <- 2 * runif(r)
opt1 <- nlm(f = log_amise_stable_log_h, p = log(h0))
opt2 <- nlm(f = log_amise_stable, p = exp(opt1$estimate))
opt3 <- nlm(f = log_amise_stable, p = h0)
expect_true(
(abs(opt1$minimum - opt2$minimum) < 1e-9) ||
(max(abs(exp(opt1$est) - abs(opt2$est))) < 1e-9)
)
expect_lt(opt1$minimum - 1e-9, opt3$minimum)
})
test_that("bw_rot_polysph() with bw0 vector and bw0 matrix", {
bw0_vec_1 <- h
bw0_vec_2 <- 0.1 * h
bw0_vec_3 <- 5 * h
bw0_mat <- rbind(bw0_vec_1, bw0_vec_2, bw0_vec_3)
expect_true(bw_rot_polysph(X = X, d = d, bw0 = bw0_vec_1,
iterlim = 1)$opt$minimum + 0.1 >=
bw_rot_polysph(X = X, d = d, bw0 = bw0_mat,
iterlim = 1)$opt$minimum)
expect_true(bw_rot_polysph(X = X, d = d, bw0 = bw0_vec_2,
iterlim = 1)$opt$minimum + 0.1 >=
bw_rot_polysph(X = X, d = d, bw0 = bw0_mat,
iterlim = 1)$opt$minimum)
expect_true(bw_rot_polysph(X = X, d = d, bw0 = bw0_vec_3,
iterlim = 1)$opt$minimum + 0.1 >=
bw_rot_polysph(X = X, d = d, bw0 = bw0_mat,
iterlim = 1)$opt$minimum)
})
## Plug-in bandwidth selectors
test_that("Same result for vMF kernel with kernel_type = 1,2", {
# Setup
r <- rpois(1, lambda = 5) + 1
d <- rpois(r, lambda = 2) + 1
h <- runif(r)
n <- 50
kappa <- 1:r
mu <- r_unif_polysph(n = 5, d = d)
X <- r_kde_polysph(n = n, X = mu, d = d, h = h)
expect_equal(
bw_rot_polysph(X = X, d = d, kernel = 1, kernel_type = 1)$bw,
bw_rot_polysph(X = X, d = d, kernel = 1, kernel_type = 2)$bw
)
})
test_that("Same result for Epa and sfp kernel with kernel_type = 1,2
and r = 1", {
# Setup
d <- rpois(1, lambda = 2) + 1
mu <- r_unif_polysph(n = 1, d = d)
X <- r_kde_polysph(n = 50, X = mu, d = d, h = runif(1))
expect_equal(
bw_rot_polysph(X = X, d = d, kernel = 2, kernel_type = 1)$bw,
bw_rot_polysph(X = X, d = d, kernel = 2, kernel_type = 2)$bw
)
expect_equal(
bw_rot_polysph(X = X, d = d, kernel = 3, kernel_type = 1)$bw,
bw_rot_polysph(X = X, d = d, kernel = 3, kernel_type = 2)$bw
)
})
test_that("Same result with kappa precomputed", {
# Setup
r <- 2
d <- rpois(r, lambda = 2) + 1
mu <- r_unif_polysph(n = 1, d = d)
X <- r_kde_polysph(n = 50, X = mu, d = d, h = runif(r))
ind <- cumsum(c(1, d + 1))
kappa <- sapply(seq_along(d), function(j) {
# Prepare data + fit vMF
data <- X[, ind[j]:(ind[j + 1] - 1)]
min(DirStats::norm2(movMF::movMF(x = data, k = 1, type = "S",
maxit = 300)$theta),
5e4) # Breaking point for later Bessels
})
expect_equal(
bw_rot_polysph(X = X + 1, d = d, kappa = kappa)$bw,
bw_rot_polysph(X = X, d = d, kappa = NULL)$bw
)
expect_equal(
bw_mrot_polysph(X = X + 1, d = d, kappa = kappa),
bw_mrot_polysph(X = X, d = d, kappa = NULL)
)
})
Try the polykde package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.