tests/testthat/test_is.R In bssm: Bayesian Inference of Non-Linear and Non-Gaussian State Space Models

context("Test importance_sample")

test_that("Test that bssm recovers the parameters of the Seatbelts model", {

#' @srrstats {G5.6, G5.6a, G5.6b, G5.7, G5.9b} Replicate Durbin&Koopman (1997)

model <- bsm_ng(Seatbelts[, "VanKilled"], distribution = "poisson",
sd_level = 1, sd_seasonal = 1, xreg = Seatbelts[, "law"],
beta = normal(0, 0, 1))

obj <- function(theta) {
model$beta[1] <- theta[1] model$R[1, 1, 1] <- theta[2]
model$R[2, 2, 1] <- theta[3] -logLik(model, particles = 0) } fit <- optim(c(0, 0, 0), obj, method = "L-BFGS-B", lower = c(-Inf, 0, 0), upper = c(10, 10, 10)) DK1997 <- c(-0.278, 0.0245, 0) # From Durbin and Koopman (1997) expect_equal(fit$par, DK1997, tol = 0.01)

# fixed seed for smooth likelihood optimization (enough only for "spdk")
fixed_seed <- sample(1:1e6, size = 1)
# Same but with importance sampling
obj <- function(theta) {
model$beta[1] <- theta[1] model$R[1, 1, 1] <- theta[2]
model$R[2, 2, 1] <- theta[3] -logLik(model, particles = 10, method = "spdk", seed = fixed_seed) } fit_is <- optim(c(0, 0, 0), obj, method = "L-BFGS-B", lower = c(-Inf, 0, 0), upper = c(10, 10, 10)) # essentially identical results in this case expect_equal(fit_is$par, DK1997, tol = 0.01)
})

test_that("Test that poisson bsm_ng give identical results with ssm_ung", {

expect_error(model_ssm_ung <- ssm_ung(y = 1:10, Z = matrix(c(1, 0), 2, 1),
T = array(c(1, 0, 1, 1), c(2, 2, 1)), R = array(diag(2, 2), c(2, 2, 1)),
a1 = matrix(0, 2, 1), P1 = diag(2, 2), state_names = c("level", "slope"),
distribution = "poisson"), NA)
expect_error(sim_ssm_ung <- importance_sample(model_ssm_ung, 4, seed = 2), NA)
expect_error(model_bsm_ng <- bsm_ng(1:10, sd_level = 2, sd_slope = 2,
P1 = diag(2, 2), distribution = "poisson"), NA)
expect_error(sim_bsm_ng <- importance_sample(model_bsm_ng, 4, seed = 2), NA)
expect_equal(sim_bsm_ng, sim_ssm_ung)
})

test_that("Test that svm still works", {
data("exchange")
model <- svm(exchange, rho = uniform(0.98, -0.999, 0.999),
sd_ar = halfnormal(0.2, 5), sigma = halfnormal(1, 2))

expect_error(sim <- importance_sample(model, 10, seed = 2), NA)

expect_gte(min(sim$weights), 0) expect_lt(max(sim$weights), Inf)
expect_true(is.finite(sum(sim$states))) expect_true(is.finite(sum(sim$weights)))
})


Try the bssm package in your browser

Any scripts or data that you put into this service are public.

bssm documentation built on May 4, 2022, 1:06 a.m.