This function runs the particle filter. The filter samples possible states (typically locations, termed particles) of an animal at each time point given the data up to (and including) that time point and a movement model.
pf_filter(
.map,
.timeline,
.state = "StateXY",
.xinit = NULL,
.xinit_pars = list(),
.yobs = list(),
.model_obs,
.model_move = move_xy(),
.n_move = 100000L,
.n_particle = 1000L,
.n_resample = as.numeric(.n_particle),
.n_record = 1000L,
.direction = c("forward", "backward"),
.verbose = getOption("patter.verbose")
)Arguments
- .map
A
SpatRasterthat defines the study area for the simulation (seeglossary). Here,.mapis used to:Simulate initial states if
.xinit = NULL(viasim_states_init());
- .timeline
A
POSIXctvector of regularly spaced time stamps that defines the timeline for the simulation. Here,.timelineis used to:Define the time steps of the simulation;
- .state, .xinit, .xinit_pars
Arguments used to simulate initial states (see
sim_states_init())..state---Acharacterthat defines theStatesub-type;.xinit---NULLor adata.tablethat defines the initial states for the simulation;.xinit_pars---A namedlistof parameters passed to the.parsargument ofsim_states_init();
- .yobs, .model_obs
Observations and observation models.
- .model_move, .n_move
The movement model.
.model_move---Acharacterstring that defines the movement model (seeModelMove);.n_move---Anintegerthat defines the number of attempts to find a legal move;
- .n_particle
An
integerthat defines the number of particles.- .n_resample
A
doublethat defines the effective sample size at which to re-sample particles.- .n_record
An
integerthat defines the number of recorded particles at each time step.- .direction
A
characterstring that defines the direction of the filter:"forward"runs the filter from.timeline[1]:.timeline[length(.timeline)];"backward"runs the filter from.timeline[length(.timeline)]:.timeline[1];
- .verbose
User output control (see
patter-progressfor supported options).
Value
The function returns a pf_particles object.
Overview
The particle filter iterates over time steps, simulating states (termed 'particles') that are consistent with the preceding data and a movement model at each time step.
The initial states for the algorithm are defined by .xinit or simulated via sim_states_init(). The word 'state' typically means location but may include additional parameters. If the initial state of the animal is known, it should be supplied via .xinit. Otherwise, sim_states_init() samples .n_particle initial coordinates from .map (via coords_init()), which are then translated into a data.table of states (that is, .xinit, via states_init()). The regions on .map from which initial coordinates are sampled can be restricted by the observations (.yobs) and other parameters. For automated handling of custom states and observation models at this stage, see the Details for sim_states_init().
The filter comprises three stages:
A movement step, in which we simulate possible states (particles) for the individual (for time steps 2, ..., T).
A weights step, in which we calculate particle weights from the log-likelihood of the data at each particle.
A re-sampling step, in which we optionally re-sample valid states using the weights.
The time complexity of the algorithm is \(O(TN)\).
The filter is implemented by the Julia function Patter.particle_filter(). To multi-thread movement and likelihood evaluations, set the number of threads via julia_connect(). See Patter.particle_filter() or JuliaCall::julia_help("particle_filter") for further information.
Algorithms
This is highly flexible routine for the reconstruction of the possible locations of an individual through time, given the data up to that time point. By modifying the observation models, it is straightforward to implement the ACPF, DCPF and ACDCPF algorithms introduced by Lavender et al. (2023) for reconstructing animal movements in passive acoustic telemetry systems using (a) acoustic time series, (b) archival time series and (c) acoustic and archival time series. pf_filter() thus replaces (and enhances) the flapper::ac(), flapper::dc(), flapper::acdc() and flapper::pf() functions.
See also
Particle filters and smoothers sample states (particles) that represent the possible locations of an individual through time, accounting for all data and the individual's movement.
To simulate artificial datasets, see
sim_*()functions (especiallysim_path_walk(),sim_array()andsim_observations()).To assemble real-world datasets for the filter, see
assemble_*()functions.pf_filter()runs the filter:To run particle smoothing, use
pf_smoother_two_filter().To map emergent patterns of space use, use a
map_*()function (such asmap_pou(),map_dens()andmap_hr()).
Examples
if (julia_run()) {
library(data.table)
library(dtplyr)
library(dplyr, warn.conflicts = FALSE)
#### Julia set up
julia_connect()
set_seed()
#### Define study period
# The resolution of the timeline defines the resolution of the simulation
timeline <- seq(as.POSIXct("2016-01-01", tz = "UTC"),
length.out = 100L, by = "2 mins")
#### Define study area
# `map` is a SpatRaster that defines the region within which movements are permitted
# Here, we consider the movements of an aquatic animal in Scotland
# ... `map` represents the bathymetry in the relevant region
# Use `set_map()` to export the map to `Julia`
map <- dat_gebco()
terra::plot(map)
set_map(map)
#### --------------------------------------------------
#### Simulation-based workflow
#### Scenario
# We have studied the movements of flapper skate off the west coast of Scotland.
# We have collected acoustic detections at receivers and depth time series.
# Here, we simulate movements and acoustic/archival observations arising from movements.
# We then apply particle filtering to the 'observations' to reconstruct
# ... simulated movements.
#### Simulate an acoustic array
moorings <- sim_array(.map = map,
.timeline = timeline,
.n_receiver = 100L)
# `moorings` includes the following default observation model parameters
# * These describe how detection probability declines with distance from a receiver
# * These are the parameters of the `ModelObsAcousticLogisTrunc` observation model
a <- moorings$receiver_alpha[1]
b <- moorings$receiver_beta[1]
g <- moorings$receiver_gamma[1]
d <- seq(1, 1000, by = 1)
plot(d, ifelse(d <= g, plogis(a * b * d), 0),
ylab = "Detection probability",
xlab = "Distance (m)",
type = "l")
#### Simulate a movement path
# Define `State` sub-type
# > We will consider the animal's movement in two-dimensions (x, y)
state <- "StateXY"
# Define the maximum moveable distance (m) between two time steps
mobility <- 750
# Define the movement model
# > We consider a two-dimensional random walk
move <-
move_xy(dbn_length = glue::glue("truncated(Gamma(1, 250.0), upper = {mobility})"),
dbn_angle = "Uniform(-pi, pi)")
# Simulate a path
paths <- sim_path_walk(.map = map,
.timeline = timeline,
.state = state,
.model_move = move)
#### Simulate observations
# Define observation model(s)
# * We simulate acoustic observations and depth time series
# * Acoustic observations are simulated according to `ModelObsAcousticLogisTrunc`
# * Depth observations are simulated according to `ModelObsDepthUniform`
models <- c("ModelObsAcousticLogisTrunc", "ModelObsDepthUniform")
# Define a `list` of parameters for the observation models
# (See `?ModelObs` for details)
pars_1 <-
moorings |>
select(sensor_id = "receiver_id", "receiver_x", "receiver_y",
"receiver_alpha", "receiver_beta", "receiver_gamma") |>
as.data.table()
pars_2 <- data.table(sensor_id = 1L,
depth_shallow_eps = 10,
depth_deep_eps = 10)
pars <- list(pars_1, pars_2)
# Simulate observational datasets
obs <- sim_observations(.timeline = timeline,
.model_obs = models,
.model_obs_pars = pars)
summary(obs)
head(obs$ModelObsAcousticLogisTrunc[[1]])
head(obs$ModelObsDepthUniform[[1]])
# Identify detections
detections <-
obs$ModelObsAcousticLogisTrunc[[1]] |>
filter(obs == 1L) |>
as.data.table()
# Collate datasets for filter
# > This requires a list of datasets, one for each data type
yobs <- list(obs$ModelObsAcousticLogisTrunc[[1]], obs$ModelObsDepthUniform[[1]])
#### Example (1): Run the filter using the default options
fwd <- pf_filter(.map = map,
.timeline = timeline,
.state = state,
.xinit = NULL, .xinit_pars = list(mobility = mobility),
.yobs = yobs,
.model_obs = models,
.model_move = move,
.n_particle = 1e4L)
## Output object structure:
# The function returns a `pf_particles` list:
class(fwd)
summary(fwd)
# `xinit` is a `data.table` of initial particles
fwd$xinit
# `states` is a `data.table` of particles:
fwd$states
# `diagnostics` is a `data.table` of filter diagnostics
fwd$diagnostics
# fwd$convergence is a Boolean that defines whether or not the algorithm converged
fwd$convergence
## Map output states:
# Map particle coordinates
pf_plot_xy(.map = map, .coord = fwd$states, .steps = 1L)
# Map a utilisation distribution
map_dens(.map = map,
.coord = fwd$states,
sigma = spatstat.explore::bw.diggle)
## Analyse filter diagnostics
# `maxlp` is the maximum log-posterior at each time step
# > exp(fwd$diagnostics$maxlp) = the highest likelihood score at each time step (0-1)
# > This should not be too low!
plot(fwd$diagnostics$timestamp, fwd$diagnostics$maxlp, type = "l")
# `ess` is the effective sample size
# > For a 2D distribution, >= 500 particles at each time step should be sufficient
# > But we often have a low ESS at the moment of a detection
# > If this is too low, we can trying boosting `.n_particle`
plot(fwd$diagnostics$timestamp, fwd$diagnostics$ess, type = "l")
points(detections$timestamp, rep(0, nrow(detections)),
pch = 21, col = "red", bg = "red", cex = 0.5)
abline(h = 500, col = "royalblue", lty = 3)
#### Example (2): Customise initial states for the filter
# Mark the known starting location on `.map`
# > Initial states are automatically sampled from `.map`
origin <- terra::setValues(map, NA)
cell <- terra::cellFromXY(map, cbind(paths$x[1], paths$y[1]))
origin[cell] <- paths$map_value[1]
fwd <- pf_filter(.map = origin,
.timeline = timeline,
.state = state,
.xinit = NULL, .xinit_pars = list(mobility = mobility),
.yobs = yobs,
.model_obs = models,
.model_move = move,
.n_particle = 1e4L)
# Specify `.xinit` manually
# > `.xinit` is resampled, as required, to generate `.n_particle` initial states
xinit <- data.table(map_value = paths$map_value[1], x = paths$x[1], y = paths$y[1])
fwd <- pf_filter(.map = map,
.timeline = timeline,
.state = state,
.xinit = xinit,
.yobs = yobs,
.model_obs = models,
.model_move = move,
.n_particle = 1e4L)
#### Example (3): Customise selected settings
# Boost the number of particles
fwd <- pf_filter(.map = map,
.timeline = timeline,
.state = state,
.xinit = xinit,
.yobs = yobs,
.model_obs = models,
.model_move = move,
.n_particle = 1e5L)
# Change the threshold ESS for resampling
fwd <- pf_filter(.map = map,
.timeline = timeline,
.state = state,
.xinit = xinit,
.yobs = yobs,
.model_obs = models,
.model_move = move,
.n_particle = 1e5L,
.n_resample = 1000L)
# Change the number of particles retained in memory
fwd <- pf_filter(.map = map,
.timeline = timeline,
.state = state,
.xinit = xinit,
.yobs = yobs,
.model_obs = models,
.model_move = move,
.n_particle = 1e5L,
.n_record = 2000L)
#### Example (2): Run the filter backwards
# (A forward and backward run is required for the two-filter smoother)
fwd <- pf_filter(.map = map,
.timeline = timeline,
.state = state,
.xinit = NULL, .xinit_pars = list(mobility = mobility),
.yobs = yobs,
.model_obs = models,
.model_move = move,
.n_particle = 1e4L,
.direction = "backward")
#### --------------------------------------------------
#### Real-world examples
#### Assemble datasets
# Define datasets for a selected animal
# > Here, we consider detections and depth time series from an example flapper skate
individual_id <- NULL
acc <- dat_acoustics[individual_id == 25L, ][, individual_id := NULL]
arc <- dat_archival[individual_id == 25L, ][, individual_id := NULL]
# Define a timeline
# * We can do this manually or use the observations to define a timeline:
timeline <- assemble_timeline(list(acc, arc), .step = "2 mins", .trim = TRUE)
timeline <- timeline[1:1440]
range(timeline)
# Assemble a timeline of acoustic observations (0, 1) and model parameters
# * The default acoustic observation model parameters are taken from `.moorings`
model_1 <- "ModelObsAcousticLogisTrunc"
acoustics <- assemble_acoustics(.timeline = timeline,
.acoustics = acc,
.moorings = dat_moorings)
# Assemble a timeline of archival observations and model parameters
# * Here, we include model parameters for `ModelObsDepthNormalTrunc`
model_2 <- "ModelObsDepthNormalTrunc"
archival <- assemble_archival(.timeline = timeline,
.archival =
arc |>
rename(obs = depth) |>
mutate(depth_sigma = 50,
depth_deep_eps = 20))
#### Visualise realisations of the movement model (the prior)
# We will use the same movement model as in previous examples
sim_path_walk(.map = map,
.timeline = timeline,
.state = state,
.model_move = move,
.n_path = 10L, .one_page = TRUE)
#### Example (1): Run filter forwards
fwd <- pf_filter(.map = map,
.timeline = timeline,
.state = state,
.xinit = NULL, .xinit_pars = list(mobility = mobility),
.yobs = list(acoustics, archival),
.model_obs = c(model_1, model_2),
.model_move = move)
#### Example (2): Run the filter backwards
bwd <- pf_filter(.map = map,
.timeline = timeline,
.state = state,
.xinit = NULL, .xinit_pars = list(mobility = mobility),
.yobs = list(acoustics, archival),
.model_obs = c(model_1, model_2),
.model_move = move,
.direction = "backward")
}
#> `patter::julia_connect()` called @ 2024-08-01 11:42:26...
#> ... Running `Julia` setup via `JuliaCall::julia_setup()`...
#> ... Validating Julia installation...
#> ... Setting up Julia project...
#> ... Handling dependencies...
#> Warning: `JULIA_NUM_THREADS` could not be set via `.threads`.
#> ... `Julia` set up with 8 thread(s).
#> `patter::julia_connect()` call ended @ 2024-08-01 11:42:31 (duration: ~5 sec(s)).
#> `patter::pf_filter()` called @ 2024-08-01 11:42:32...
#> ... 11:42:32: Checking user inputs...
#> ... 11:42:32: Setting initial states...
#> ... 11:42:32: Setting observations...
#> ... 11:42:32: Running filter...
#> ... 11:42:33: Collating outputs...
#> `patter::pf_filter()` call ended @ 2024-08-01 11:42:33 (duration: ~1 sec(s)).
#> `patter::map_dens()` called @ 2024-08-01 11:42:33...
#> ... 11:42:33: Processing `.map`...
#> ... 11:42:33: Building XYM...
#> ... 11:42:34: Defining `ppp` object...
#> Observation window is gridded.
#> ... 11:42:34: Estimating density surface...
#> ... 11:42:43: Scaling density surface...
#> `patter::map_dens()` call ended @ 2024-08-01 11:42:43 (duration: ~10 sec(s)).
#> `patter::pf_filter()` called @ 2024-08-01 11:42:43...
#> ... 11:42:43: Checking user inputs...
#> ... 11:42:43: Setting initial states...
#> Warning: `terra::spatSample()` failed: implementing exhaustive sampling...
#> Warning: [spatSample] fewer samples than requested are available
#> ... 11:42:43: Setting observations...
#> ... 11:42:43: Running filter...
#> ... 11:42:43: Collating outputs...
#> `patter::pf_filter()` call ended @ 2024-08-01 11:42:43 (duration: ~0 sec(s)).
#> `patter::pf_filter()` called @ 2024-08-01 11:42:43...
#> ... 11:42:43: Checking user inputs...
#> ... 11:42:43: Setting initial states...
#> ... 11:42:43: Setting observations...
#> ... 11:42:43: Running filter...
#> ... 11:42:44: Collating outputs...
#> `patter::pf_filter()` call ended @ 2024-08-01 11:42:44 (duration: ~1 sec(s)).
#> `patter::pf_filter()` called @ 2024-08-01 11:42:44...
#> ... 11:42:44: Checking user inputs...
#> ... 11:42:44: Setting initial states...
#> ... 11:42:44: Setting observations...
#> ... 11:42:44: Running filter...
#> ... 11:42:47: Collating outputs...
#> `patter::pf_filter()` call ended @ 2024-08-01 11:42:47 (duration: ~3 sec(s)).
#> `patter::pf_filter()` called @ 2024-08-01 11:42:47...
#> ... 11:42:47: Checking user inputs...
#> ... 11:42:47: Setting initial states...
#> ... 11:42:47: Setting observations...
#> ... 11:42:47: Running filter...
#> ... 11:42:48: Collating outputs...
#> `patter::pf_filter()` call ended @ 2024-08-01 11:42:48 (duration: ~1 sec(s)).
#> `patter::pf_filter()` called @ 2024-08-01 11:42:48...
#> ... 11:42:48: Checking user inputs...
#> ... 11:42:48: Setting initial states...
#> ... 11:42:48: Setting observations...
#> ... 11:42:48: Running filter...
#> ... 11:42:51: Collating outputs...
#> `patter::pf_filter()` call ended @ 2024-08-01 11:42:51 (duration: ~3 sec(s)).
#> `patter::pf_filter()` called @ 2024-08-01 11:42:51...
#> ... 11:42:51: Checking user inputs...
#> ... 11:42:51: Setting initial states...
#> ... 11:42:51: Setting observations...
#> ... 11:42:51: Running filter...
#> ... 11:42:51: Collating outputs...
#> `patter::pf_filter()` call ended @ 2024-08-01 11:42:51 (duration: ~0 sec(s)).
#> `patter::pf_filter()` called @ 2024-08-01 11:42:52...
#> ... 11:42:52: Checking user inputs...
#> ... 11:42:52: Setting initial states...
#> ... 11:42:52: Setting observations...
#> ... 11:42:53: Running filter...
#> ... 11:42:53: Collating outputs...
#> `patter::pf_filter()` call ended @ 2024-08-01 11:42:53 (duration: ~1 sec(s)).
#> `patter::pf_filter()` called @ 2024-08-01 11:42:53...
#> ... 11:42:53: Checking user inputs...
#> ... 11:42:53: Setting initial states...
#> ... 11:42:53: Setting observations...
#> ... 11:42:53: Running filter...
#> ... 11:42:54: Collating outputs...
#> `patter::pf_filter()` call ended @ 2024-08-01 11:42:54 (duration: ~1 sec(s)).