tests/testthat/test-real-pathology-fixes.R

In bayesianETAS: Bayesian Estimation of the Temporal and Spatio-Temporal ETAS Models for Earthquake Occurrences

test_that("duplicate event times are rejected by default and deterministic jitter is opt-in", {
  ts <- c(0.1, 0.1, 0.35, 0.8, 1.4)
  ms <- c(3.0, 3.1, 3.2, 3.1, 3.4)
  init <- c(0.6, 0.2, 0.4, 0.08, 1.5)

  expect_error(
    maxLikelihoodETAS(ts, ms, 3, maxTime = 2, initval = init),
    "duplicate event times detected"
  )

  fit <- maxLikelihoodETAS(ts, ms, 3, maxTime = 2, initval = init,
                           handle_ties = "jitter")
  expect_true(isTRUE(attr(fit, "tie_handling")$jittered))
  expect_gt(fit$params["c"], attr(fit, "c_lower"))
  expect_true(is.finite(fit$loglik))
})

test_that("c_lower is enforced in MLE and MCMC", {
  ts <- c(0.1, 0.3, 0.7, 1.2, 1.8, 2.4)
  ms <- c(3.0, 3.1, 3.2, 3.1, 3.4, 3.3)
  bad_init <- c(0.5, 0.2, 0.4, 0.05, 1.5)
  good_init <- c(0.5, 0.2, 0.4, 0.2, 1.5)

  expect_error(
    maxLikelihoodETAS(ts, ms, 3, maxTime = 3, initval = bad_init, c_lower = 0.1),
    "c must be greater than c_lower"
  )

  fit <- maxLikelihoodETAS(ts, ms, 3, maxTime = 3,
                           initval = good_init, c_lower = 0.1)
  expect_gt(fit$params["c"], 0.1)

  set.seed(42)
  samples <- estimateETAS(ts, ms, 3, maxTime = 3, sims = 8, burnin = 0,
                          initval = good_init, c_lower = 0.1)
  expect_true(all(samples$c > 0.1))
})

test_that("reported MLE log-likelihood equals likelihood at returned parameters", {
  ts <- c(0.1, 0.25, 0.55, 0.9, 1.4, 2.0)
  ms <- c(3.0, 3.2, 3.1, 3.4, 3.3, 3.5)
  fit <- maxLikelihoodETAS(ts, ms, 3, maxTime = 2.5,
                           initval = c(0.6, 0.2, 0.4, 0.08, 1.5))
  ll <- bayesianETAS:::ETASlikelihood(
    ts, ms, 2.5, 3,
    fit$params["mu"], fit$params["K"], fit$params["alpha"],
    fit$params["c"], fit$params["p"],
    c_lower = attr(fit, "c_lower")
  )

  expect_equal(as.numeric(fit$loglik), as.numeric(ll), tolerance = 1e-7)
})

test_that("c and p move in a small temporal MCMC run", {
  ts <- c(0.1, 0.22, 0.35, 0.6, 0.95, 1.4, 2.1, 2.8)
  ms <- c(3.0, 3.3, 3.1, 3.4, 3.2, 3.5, 3.1, 3.3)

  set.seed(123)
  samples <- estimateETAS(ts, ms, 3, maxTime = 3.5, sims = 20, burnin = 0,
                          initval = c(0.8, 0.2, 0.4, 0.08, 1.5),
                          c_lower = 0.005)

  expect_gt(length(unique(round(samples$c, 10))), 1)
  expect_gt(length(unique(round(samples$p, 10))), 1)
  expect_true(all(is.finite(as.matrix(samples))))
  expect_true(all(samples$c > 0.005))
  expect_true(all(samples$p > 1))
})

test_that("synthetic temporal and spatio-temporal fits return finite legal parameters", {
  ts <- c(0.1, 0.25, 0.5, 0.85, 1.3, 1.9, 2.7)
  ms <- c(3.0, 3.2, 3.1, 3.4, 3.3, 3.5, 3.2)
  xs <- c(0.0, 0.2, 0.35, 0.6, 0.8, 1.0, 1.3)
  ys <- c(0.0, 0.1, 0.25, 0.35, 0.6, 0.7, 0.95)

  tfit <- maxLikelihoodETAS(ts, ms, 3, maxTime = 3,
                            initval = c(0.8, 0.2, 0.4, 0.08, 1.5),
                            c_lower = 0.005)
  expect_true(all(is.finite(tfit$params)))
  expect_gte(tfit$params["mu"], 0)
  expect_gte(tfit$params["K"], 0)
  expect_gte(tfit$params["alpha"], 0)
  expect_gt(tfit$params["c"], 0.005)
  expect_gt(tfit$params["p"], 1)

  sfit <- maxLikelihoodETAS(ts, ms, 3, maxTime = 3, xs = xs, ys = ys,
                            area = 2, c_lower = 0.005,
                            initval = c(0.8, 0.2, 0.4, 0.08, 1.5, 0.3, 0.3))
  expect_true(all(is.finite(sfit$params)))
  expect_gte(sfit$params["mu"], 0)
  expect_gte(sfit$params["K"], 0)
  expect_gte(sfit$params["alpha"], 0)
  expect_gt(sfit$params["c"], 0.005)
  expect_gt(sfit$params["p"], 1)
  expect_gt(sfit$params["sigma2x"], 0)
  expect_gt(sfit$params["sigma2y"], 0)
})