tests/testthat/test-compose-sum-ar1-matern.R

In ngme2: Linear Latent Non-Gaussian Models with Flexible Distributions

test_that("test sum of ar1 and matern", {
  skip_on_cran()
  withr::local_seed(123)
  library(ngme2)
  library(fmesher)

  # 1. Simulation
  n_obs <- 1000
  time_points <- 1:n_obs

  # AR1 process
  rho_true <- 0.7
  ar1_model <- f(time_points, model = ar1(rho = rho_true), noise = noise_normal(sigma = 1.5))
  W_ar1 <- simulate(ar1_model)[[1]]
  sd(W_ar1)

  # Matern process
  # Create a 1D mesh
  mesh <- fm_mesh_1d(seq(1, n_obs, length.out = 50))
  kappa_true <- 0.3
  matern_model <- f(time_points, model = matern(mesh, kappa = kappa_true), noise = noise_normal(sigma = 4.0))
  W_matern <- simulate(matern_model)[[1]]
  sd(W_matern)

  # Measurement noise
  sigma_eps_true <- 0.5
  Y <- W_ar1 + W_matern + rnorm(n_obs, sd = sigma_eps_true)

  # 2. Estimation
  # We need to specify the model structure for estimation
  # Note: For Matern 1D, we need to provide the mesh

  out <- ngme(
    Y ~ 0 + f(time_points, model = ar1(), name = "ar1_comp") +
      f(time_points, model = matern(mesh), name = "matern_comp"),
    data = data.frame(Y = Y),
    control_opt = control_opt(
      seed = 123,
      iterations = 1000,
      optimizer = adamW(),
      n_parallel_chain = 4,
      print_check_info = FALSE,
      verbose = FALSE
    )
  )
  traceplot(out, name = "ar1_comp", hline = c(rho_true, 1.5))
  traceplot(out, name = "matern_comp", hline = c(kappa_true, 4))

  # 3. Comparison
  res_ar1 <- ngme_result(out, "ar1_comp")
  res_matern <- ngme_result(out, "matern_comp")
  res_noise <- ngme_result(out, "data") # Measurement noise is usually under "data" or similar if not specified in f()

  # Parallel stochastic fitting is somewhat noisy here, so keep the assertions
  # focused on recovering a sensible decomposition rather than exact targets.
  expect_true(is.finite(res_ar1$rho))
  expect_true(res_ar1$rho > 0.4 && res_ar1$rho < 1.05)
  expect_true(is.finite(res_matern$kappa))
  expect_true(res_matern$kappa > 0)
  expect_true(is.finite(res_noise$sigma))
  expect_true(res_noise$sigma > 0.2 && res_noise$sigma < 0.85)

  # Print results for manual inspection
  print(paste("True AR1 rho:", rho_true, "Est:", res_ar1$rho))
  print(paste("True Matern kappa:", kappa_true, "Est:", res_matern$kappa)) # Assuming kappa is directly available or we need to extract it
  print(paste("True Sigma eps:", sigma_eps_true, "Est:", res_noise$sigma))
})