tests/testthat/test_SFPCA.R

In DistributionalDDA: Distributional Data Analysis in Bayes Spaces

# Regression tests for the SFPCA loading back-transform (fixed 2026-07-02):
# with S = M t(M) (M = t(chol(S))) and whitened coefficients C M, a whitened
# eigenvector v_tilde pulls back to the loading v = t(M)^{-1} v_tilde.  The
# implementation previously used M^{-1} v_tilde, which is wrong for the
# non-symmetric Cholesky factor: eigenfunctions were not S-orthonormal and
# the modes-of-variation displays were rotated.

test_that("SFPCA loadings are S-orthonormal and consistent with scores", {
  set.seed(7)
  basis <- B2Unif01Basis7
  S <- getMetricTensor(basis)
  x <- BDataSplines(matrix(rnorm(30 * 7, sd = 0.3), nrow = 30), basis = basis)
  p <- suppressMessages(princomp(x))
  V <- unclass(p$loadings)
  # S-orthonormality of the back-transformed loadings
  expect_equal(t(V) %*% S %*% V, diag(ncol(V)), tolerance = 1e-10,
               ignore_attr = TRUE)
  # scores equal the S-inner products of centered data with the loadings
  Cc <- sweep(unclass(x), 2, p$center, "-")
  expect_equal(unname(Cc %*% S %*% V), unname(p$scores), tolerance = 1e-10)
  # predict() on the training data reproduces the scores
  expect_equal(unname(predict(p, x)), unname(p$scores), tolerance = 1e-10)
})

test_that("multivariate SFPCA loadings are block-S-orthonormal", {
  set.seed(8)
  basis <- B2Unif01Basis5
  S <- getMetricTensor(basis)
  x <- BDataSplines(matrix(rnorm(40 * 5, sd = 0.3), nrow = 40), basis = basis)
  g <- rep(c("a", "b"), each = 20)
  pm <- suppressMessages(princomp(x, group = g, multivariate = TRUE))
  V <- unclass(pm$loadings)
  Sstack <- as.matrix(Matrix::bdiag(S, S))
  expect_equal(t(V) %*% Sstack %*% V, diag(ncol(V)), tolerance = 1e-10,
               ignore_attr = TRUE)
})