test_that("Hawkes simulation (1)", {
## sim_hawkes example
set.seed(1234)
times <- sim_hawkes(mu = 0.3, alpha = 4, beta = 5)
expect_equal(round(times[1:6], 2),
c(4.15, 4.95, 7.80, 14.64, 18.69, 33.61),
tolerance = 0.01)
})
test_that("Hawkes simulation (2)", {
## sim_hawkes example
set.seed(1234)
times <- sim_hawkes(mu = 0.3, alpha = 4, beta = 5, method = "2")
expect_equal(round(times[1:6], 2),
c(4.15, 4.22, 4.40, 4.42, 4.43, 5.54),
tolerance = 0.01)
})
test_that("Simple Hawkes model fitting", {
## NIWA retweets
data(retweets_niwa, package = "stelfi")
times <- unique(sort(as.numeric(difftime(retweets_niwa,
min(retweets_niwa),
units = "mins"))))
params <- c(mu = 9, alpha = 3, beta = 10)
fit <- fit_hawkes(times = times, parameters = params)
pars <- as.numeric(get_coefs(fit)[, 1])
expect_equal(pars[1], 0.06328, tolerance = 1)
expect_equal(pars[2], 0.07597, tolerance = 1)
expect_equal(pars[3], 0.07911, tolerance = 1)
})
test_that("Non-homogeneous Hawkes model fitting", {
set.seed(1)
times <- sim_hawkes(mu = 0.3, alpha = 4, beta = 5)
background <- function(params, times) {
A <- exp(params[[1]])
B <- stats::plogis(params[[2]]) * A
return(A + B * sin(times))
}
background_integral <- function(params, x) {
A <- exp(params[[1]])
B <- stats::plogis(params[[2]]) * A
return((A * x) - B * cos(x))
}
param <- list(alpha = 0.5, beta = 1.5)
background_param <- list(1, 1)
fit <- fit_hawkes_cbf(times = times, parameters = param,
background = background,
background_integral = background_integral,
background_parameters = background_param)
estA <- exp(fit$background_parameters[1])
estB <- plogis(fit$background_parameters[2]) * exp(fit$background_parameters[1])
expect_equal(estA, 0.2387370, tolerance = 0.1)
expect_equal(estB, 0.02600446, tolerance = 0.1)
})
test_that("Multivariate Hawkes model fitting", {
data(multi_hawkes, package = "stelfi")
fit <- stelfi::fit_mhawkes(times = multi_hawkes$times, stream = multi_hawkes$stream,
parameters = list(mu = c(0.2,0.2),
alpha = matrix(c(0.5,0.1,0.1,0.5),byrow = TRUE,nrow = 2),
beta = c(0.7,0.7)))
ll <- fit$fn()
expect_equal(ll, 132, tolerance = 2)
})
test_that("LGCP model fitting (spatial)", {
if(requireNamespace("fmesher")){
data(xyt, package = "stelfi")
domain <- sf::st_as_sf(xyt$window)
locs <- data.frame(x = xyt$x, y = xyt$y)
bnd <- fmesher::fm_as_segm(as.matrix(sf::st_coordinates(domain)[, 1:2]))
smesh <- fmesher::fm_mesh_2d(boundary = bnd,
max.edge = 0.75, cutoff = 0.3)
fit <- fit_lgcp(locs = locs, sf = domain, smesh = smesh,
parameters = c(beta = 0, log_tau = log(1),
log_kappa = log(1)))
pars <- as.numeric(get_coefs(fit)[, 1])
expect_equal(pars[1], 2.45, tolerance = 0.5)
expect_equal(pars[2], -1.32, tolerance = 0.1)
expect_equal(pars[3], 0.95, tolerance = 0.1)
}
})
test_that("Simulate LGCP (spatial)", {
if(requireNamespace("fmesher")){
data(xyt, package = "stelfi")
domain <- sf::st_as_sf(xyt$window)
bnd <- fmesher::fm_as_segm(as.matrix(sf::st_coordinates(domain)[, 1:2]))
smesh <- fmesher::fm_mesh_2d(boundary = bnd, max.edge = 0.75, cutoff = 0.3)
parameters <- c(beta = 1, log_tau = log(1), log_kappa = log(1))
set.seed(91234)
sim <- sim_lgcp(parameters = parameters, sf = domain, smesh = smesh)
expect_equal(round(c(sim$x[1:3]), 3),
c(0.172, 0.175, 0.298),
tolerance = 0.5)
expect_equal(length(sim$x),length(sim$y))
}
})
test_that("LGCP model fitting (spatiotemporal)", {
skip_on_cran()
if(requireNamespace("fmesher")){
data(xyt, package = "stelfi")
domain <- sf::st_as_sf(xyt$window)
ndays <- 2
locs <- data.frame(x = xyt$x, y = xyt$y, t = xyt$t)
bnd <- fmesher::fm_as_segm(as.matrix(sf::st_coordinates(domain)[, 1:2]))
w0 <- 2
smesh <- fmesher::fm_mesh_2d(boundary = bnd,
max.edge = 0.75, cutoff = 0.3)
tmesh <- fmesher::fm_mesh_1d(seq(0, ndays, by = w0))
fit <- fit_lgcp(locs = locs, sf = domain, smesh = smesh, tmesh = tmesh,
parameters = c(beta = 0, log_tau = log(1),
log_kappa = log(1), atanh_rho = 0.2))
pars <- as.numeric(get_coefs(fit)[, 1])
expect_equal(pars[1], 0.31, tolerance = 1)
expect_equal(pars[2], 1.68, tolerance = 1)
expect_equal(pars[3], -1.06, tolerance = 1)
}
})
test_that("Simulate LGCP (spatiotemporal)", {
skip_on_cran()
if(requireNamespace("fmesher")){
data(xyt, package = "stelfi")
domain <- sf::st_as_sf(xyt$window)
bnd <- fmesher::fm_as_segm(as.matrix(sf::st_coordinates(domain)[, 1:2]))
smesh <- fmesher::fm_mesh_2d(boundary = bnd, max.edge = 0.75, cutoff = 0.3)
set.seed(91234)
ndays <- 2
w0 <- 2
tmesh <- fmesher::fm_mesh_1d(seq(0, ndays, by = w0))
parameters <- c(beta = 1, log_tau = log(1), log_kappa = log(1), atanh_rho = 0.2)
sim <- sim_lgcp(parameters = parameters, sf = domain, smesh = smesh, tmesh = tmesh)
expect_equal(round(c(sim$x[1:3, 2]), 3), c(0.252, 0.674, 0.575), tolerance = 0.1)
expect_equal(c(sim$y[1:3]), c(0, 0, 1), tolerance = 0.1)
}
})
test_that("LGCP model fitting (marked)", {
skip_on_cran()
if(requireNamespace("fmesher")){
data(marked, package = "stelfi")
loc.d <- 3 * cbind(c(0, 1, 1, 0, 0), c(0, 0, 1, 1, 0))
domain <- sf::st_sf(geometry = sf::st_sfc(sf::st_polygon(list(loc.d))))
smesh <- fmesher::fm_mesh_2d(loc.domain = loc.d, offset = c(0.3, 1),
max.edge = c(0.3, 0.7), cutoff = 0.05)
locs <- cbind(x = marked$x, y = marked$y)
marks <- cbind(m1 = marked$m1) ## Gaussian
parameters <- list(betamarks = matrix(0, nrow = 1, ncol = ncol(marks)),
log_tau = log(1), log_kappa = log(1),
marks_coefs_pp = rep(0, ncol(marks)), betapp = 0)
fit <- fit_mlgcp(locs = locs, marks = marks,
sf = domain, smesh = smesh,
parameters = parameters, methods = 0,
fields = 0)
pars <- as.numeric(get_coefs(fit)[, 1])
expect_equal(pars[1], 9.90, tolerance = 0.1)
expect_equal(pars[3], -0.279, tolerance = 0.1)
}
})
test_that("LGCP model fitting (marked) with covariate overlap", {
skip_on_cran()
if(requireNamespace("fmesher")){
data(marked, package = "stelfi")
loc.d <- 3 * cbind(c(0, 1, 1, 0, 0), c(0, 0, 1, 1, 0))
domain <- sf::st_sf(geometry = sf::st_sfc(sf::st_polygon(list(loc.d))))
smesh <- fmesher::fm_mesh_2d(loc.domain = loc.d, offset = c(0.3, 1),
max.edge = c(0.3, 0.7), cutoff = 0.05)
locs <- cbind(x = marked$x, y = marked$y)
marks <- cbind(m1 = marked$m1) ## Gaussian
set.seed(132)
covs <- cbind(cov = rnorm(smesh$n))
parameters <- list(betamarks = matrix(0, nrow = 2, ncol = ncol(marks)),
log_tau = log(1), log_kappa = log(1),
marks_coefs_pp = rep(0, ncol(marks)), betapp = c(0, 1))
fit <- fit_mlgcp(locs = locs, marks = marks,
sf = domain, smesh = smesh,
parameters = parameters, methods = 0,
fields = 0, covariates = covs,
pp_covariates = 1,
marks_covariates = 1)
pars <- as.numeric(get_coefs(fit)[, 1])
expect_equal(pars[4], 0.08150350, tolerance = 0.1)
}
})
test_that("Spatial self-exciting (no GMRF)", {
skip_on_cran()
if(requireNamespace("fmesher")){
data(xyt, package = "stelfi")
N <- 50
locs <- data.frame(x = xyt$x[1:N], y = xyt$y[1:N])
times <- xyt$t[1:N]
domain <- sf::st_as_sf(xyt$window)
bnd <- fmesher::fm_as_segm(as.matrix(sf::st_coordinates(domain)[, 1:2]))
smesh <- fmesher::fm_mesh_2d(boundary = bnd,
max.edge = 0.75, cutoff = 0.3)
param <- list(mu = 3, alpha = 1, beta = 3, xsigma = 0.2,
ysigma = 0.2, rho = 0.8)
fit <- fit_stelfi(times = times, locs = locs, sf = domain, smesh = smesh, parameters = param)
pars <- as.numeric(get_coefs(fit)[, 1])
expect_equal(pars[1], 0.3, tolerance = 0.2)
expect_equal(pars[2], -1.1, tolerance = 0.2)
expect_equal(pars[7], -0.17, tolerance = 0.2)
}
})
test_that("Spatial self-exciting (GMRF)", {
skip_on_cran()
if(requireNamespace("fmesher")){
data(xyt, package = "stelfi")
N <- 50
locs <- data.frame(x = xyt$x[1:N], y = xyt$y[1:N])
times <- xyt$t[1:N]
domain <- sf::st_as_sf(xyt$window)
bnd <- fmesher::fm_as_segm(as.matrix(sf::st_coordinates(domain)[, 1:2]))
smesh <- fmesher::fm_mesh_2d(boundary = bnd,
max.edge = 0.75, cutoff = 0.3)
param <- list( mu = 5, alpha = 1, beta = 3, kappa = 0.9, tau = 1, xsigma = 0.2, ysigma = 0.2, rho = 0.8)
fit <- fit_stelfi(times = times, locs = locs, sf = domain, smesh = smesh, parameters = param, GMRF = TRUE)
pars <- as.numeric(get_coefs(fit)[, 1])
expect_equal(pars[1], 0.003, tolerance = 0.02)
expect_equal(pars[2], -5.88, tolerance = 0.2)
}
})
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.