tests/testthat/test_DomainScale.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

# Regression tests for the domain-width scale bug (fixed 2026-07-02):
# BayesSpaceBasis.basisfd stored the Gram matrix as the raw Lebesgue integral
# int zeta_i zeta_j dlambda over rangeval, while every inner product in the
# package (CovarianceOperator, expectation, scalarProductOfLogDensities) is a
# MEAN with respect to the uniform reference measure. On the unit interval the
# two conventions coincide; on a domain of width W the projection
# solve(S, <f, zeta>) returned coefficients too small by the factor W, so that
# densities converted through as.BDataSplines were flattened towards the
# reference measure (e.g., by a factor 100 for age distributions on [0, 100]).
# S is now divided by diff(rangeval).

test_that("in-span log-densities are recovered on a non-unit domain", {
  set.seed(11)
  b <- BayesSpaceBasis(BayesSpace(ReferenceMeasureUniform(0, 100, sample= runif(50000, 0, 100))),
                       nbasis = 5)
  c0 <- rnorm(5)
  f  <- BDataLogDensities(list(function(x) c(basisEval(b, x) %*% c0)), basis = b)
  cHat <- unclass(as.BDataSplines(f))
  # compare reconstructed function values on a grid (robust to the mild
  # ill-conditioning of the coefficient representation)
  xg <- seq(1, 99, length.out = 200)
  fTrue <- c(basisEval(b, xg) %*% c0)
  fHat  <- c(basisEval(b, xg) %*% c(cHat))
  expect_lt(sqrt(mean((fHat - fTrue)^2)) / sqrt(mean(fTrue^2)), 0.05)
})

test_that("Gram-based and Monte-Carlo norms agree on a non-unit domain", {
  set.seed(12)
  b <- BayesSpaceBasis(BayesSpace(ReferenceMeasureUniform(0, 100, sample= runif(50000,0, 100))),
                       nbasis = 5)
  xs <- BDataSplines(matrix(rnorm(10), nrow = 2), basis = b)
  fast <- norm(xs)
  slow <- norm(as.BDataLogDensities(xs))
  expect_equal(fast, slow, tolerance = 0.05)
})

test_that("scalarProduct on splines matches the functional route off [0,1]", {
  set.seed(13)
  b <- BayesSpaceBasis(BayesSpace(ReferenceMeasureUniform(0, 50, sample= runif(5000,0, 50))),
                       nbasis = 5)
  xs <- BDataSplines(matrix(rnorm(10), nrow = 2), basis = b)
  spSpl <- scalarProduct(xs, xs)
  spFun <- scalarProduct(as.BDataLogDensities(xs), as.BDataLogDensities(xs))
  expect_equal(unname(as.matrix(spSpl)), unname(as.matrix(spFun)),
               tolerance = 0.1)
})