tests/testthat/test-phmc.R

In pHMC: Proximal Hamiltonian Monte Carlo for Non-Smooth Bayesian Inference

test_that("Proximal operators compute correct values", {
  ## L1 soft-thresholding

  x <- c(-2, -0.5, 0, 0.5, 2)
  res_l1 <- prox_l1(x, tau = 1)
  expect_equal(res_l1, c(-1, 0, 0, 0, 1))
  expect_true(all(is.finite(res_l1)))

  ## L2 shrinkage
  x2 <- c(3, 4)
  res_l2 <- prox_l2(x2, tau = 2.5)
  expect_equal(res_l2, c(1.5, 2.0))
  expect_true(all(is.finite(res_l2)))

  ## Nuclear norm thresholding
  mat <- matrix(c(3, 0, 0, 4), nrow = 2)
  res_nuc <- prox_nuclear(mat, tau = 1)
  expect_equal(unname(as.matrix(res_nuc)),
               matrix(c(2, 0, 0, 3), nrow = 2))
  expect_true(all(is.finite(res_nuc)))

  ## Elastic net
  res_en <- prox_elastic_net(c(2, -2), tau = 0.5, alpha = 0.5)
  expect_equal(length(res_en), 2L)
  expect_true(all(is.finite(res_en)))
})

test_that("Moreau-Yosida envelope gradient is accurate", {
  x <- c(2, -2)
  grad_my <- grad_my_envelope(x, prox_fn = "l1", lambda_g = 0.1)
  expect_equal(length(grad_my), 2L)
  expect_true(all(is.finite(grad_my)))
  expect_false(any(is.na(grad_my)))
  expect_false(any(is.nan(grad_my)))

  ## 'none' should give zero gradient
  grad_none <- grad_my_envelope(x, prox_fn = "none", lambda_g = 0.1)
  expect_equal(grad_none, c(0, 0))
})

test_that("phmc runs without NAs/NaNs/warnings on 1D Bayesian Lasso", {
  set.seed(123)
  y_data <- rnorm(50, mean = 1, sd = 0.5)
  f_smooth <- function(x, y) 0.5 * sum((y - x)^2)
  grad_f_smooth <- function(x, y) -sum(y - x)

  fit <- suppressWarnings(phmc(
    fn = f_smooth,
    grad_f = grad_f_smooth,
    prox_fn = "l1",
    start = 0.5,
    data = y_data,
    lambda_g = 0.01,
    n_draws = 100,
    burnin = 20,
    epsilon = 0.01,
    L = 5
  ))

  expect_s3_class(fit, "phmc")
  expect_true(fit$accept_rate > 0)
  expect_equal(nrow(fit$draws), 80L)

  ## No NAs/NaNs anywhere
  expect_false(any(is.na(fit$draws)))
  expect_false(any(is.nan(fit$draws)))
  expect_false(any(is.na(fit$estimates)))
  expect_false(any(is.nan(fit$estimates)))
  expect_true(is.finite(fit$logLik))
  expect_true(is.finite(fit$AIC))
  expect_true(is.finite(fit$BIC))
  expect_true(is.finite(fit$DIC))

  ## S3 methods work
  co <- coef(fit)
  expect_true(is.numeric(co))
  expect_equal(length(co), 1L)
  expect_true(all(is.finite(co)))

  vc <- vcov(fit)
  expect_true(is.matrix(vc))
  expect_true(all(is.finite(vc)))

  ll <- logLik(fit)
  expect_true(is.numeric(ll))
  expect_true(is.finite(ll))

  sm <- summary(fit)
  expect_s3_class(sm, "summary.phmc")
})

test_that("phmc_tune selects a valid lambda_g", {
  set.seed(123)
  f_smooth <- function(x) 0.5 * sum(x^2)
  tune_res <- phmc_tune(
    fn = f_smooth,
    prox_fn = "l1",
    start = c(1, 1),
    lambda_grid = c(0.001, 0.01, 0.1),
    epsilon = 0.001,
    L = 5
  )

  expect_s3_class(tune_res, "phmc_tune")
  expect_true(tune_res$optimal_lambda_g %in% c(0.001, 0.01, 0.1))
  expect_true(is.data.frame(tune_res$grid_results))
  expect_false(any(is.na(tune_res$grid_results$R_lambda_g)))
  expect_false(any(is.nan(tune_res$grid_results$R_lambda_g)))
})

test_that("phmc_logistic works without NAs on simulated binary data", {
  set.seed(456)
  X <- matrix(rnorm(60), nrow = 30, ncol = 2)
  y <- rbinom(30, 1, 0.5)

  fit <- suppressWarnings(phmc_logistic(X, y, alpha = 0.5,
    lambda_g = 0.01, n_draws = 50, burnin = 10))
  expect_s3_class(fit, "phmc")
  expect_equal(length(coef(fit)), 2L)
  expect_false(any(is.na(fit$draws)))
  expect_false(any(is.nan(fit$draws)))
})

test_that("phmc_lasso works without NAs on simulated data", {
  set.seed(101)
  X <- matrix(rnorm(60), nrow = 30, ncol = 2)
  y <- X %*% c(1, -0.5) + rnorm(30, sd = 0.3)

  fit <- suppressWarnings(phmc_lasso(X, y, alpha = 0.5,
    lambda_g = 0.01, n_draws = 50, burnin = 10))
  expect_s3_class(fit, "phmc")
  expect_equal(length(coef(fit)), 2L)
  expect_false(any(is.na(fit$draws)))
})

test_that("phmc_matrix works without NAs on small matrix", {
  set.seed(789)
  Y <- matrix(rnorm(9), nrow = 3, ncol = 3)
  fit <- suppressWarnings(phmc_matrix(Y, alpha = 0.5,
    lambda_g = 0.001, n_draws = 40, burnin = 10))
  expect_s3_class(fit, "phmc")
  expect_equal(nrow(fit$draws), 30L)
  expect_false(any(is.na(fit$draws)))
})

test_that("phmc with finite-difference gradient (grad_f = NULL)", {
  set.seed(222)
  y_data <- rnorm(30, mean = 2, sd = 0.5)
  f_smooth <- function(x, y) 0.5 * sum((y - x)^2)

  fit <- suppressWarnings(phmc(
    fn = f_smooth,
    grad_f = NULL,
    prox_fn = "l1",
    start = 0.0,
    data = y_data,
    lambda_g = 0.01,
    n_draws = 50,
    burnin = 10,
    epsilon = 0.005,
    L = 5
  ))

  expect_s3_class(fit, "phmc")
  expect_false(any(is.na(fit$draws)))
  expect_false(any(is.nan(fit$draws)))
  expect_true(is.finite(fit$logLik))
})