Observation models

Source: R/models.R

ModelObs is Abstract Type in Patter.jl that groups observation model sub-types.

Details

Observation model sub-types are Julia structures that hold the parameters of observation models. From an R-user perspective, you can think of a ModelObs sub-type as an S4-class-like object, with slots for the parameters of an observation model. With an observation model structure, we can simulate new observations and evaluate the log-probability of existing observations.

The following observation models are built in to Patter.jl:

  • ModelObsAcousticLogisTrunc

  • ModelObsDepthUniform

  • ModelObsDepthNormalTrunc

See Patter.jl or JuliaCall::julia_help("ModelObs") for the fields of the built-in sub-types.

In patter, observation models are required:

Observation model sub-types should be specified as a character vector alongside a list of data.table that contain the parameter values for each model. Internally, observation model sub-types and parameters are instantiated and used to simulate observations or in the particle filter. The simulation of observations is implemented via Patter.simulate_obs(). In the particle filter, log-probabilities are evaluated by Patter.logpdf_obs(). These are generic functions. Different methods are dispatched according to the input model. For the built-in ModelObs sub-types, corresponding methods for these routines are also built-in. For custom ModelObs sub-types, the methods need to be provided.

To use custom ModelObs sub-types, see Examples.

See also

The routines in patter for the simulation of individual movements, observations and statistical modelling are built upon three Abstract Types defined in Julia:

  • See State for individual-state (location) structures;

  • See ModelMove for movement model structures;

  • See ModelObs for observation model structures;

Author

Edward Lavender

Examples

# Patter contains multiple built-in `ModelObs` sub-types that you can use
# ... (with custom parameters) simulate observations and for particle filtering.
# To use a new sub-type, follow the workflow below. Some extra work is required
# ... because we have to register the sub-type in `Julia` and write the
# ... required methods to simulate observations and/or calculate log probabilities.

if (julia_run()) {

  library(JuliaCall)
  library(data.table)

  #### Julia set up
  # Connect to Julia
  julia <- julia_connect()
  # Set seed
  set_seed()
  # Export map to Julia
  map <- dat_gebco()
  set_map(map)

  #### Simulate path(s)
  # > We simulate a path in four dimensions (see `?StateXYZD`)
  timeline <- seq(as.POSIXct("2016-01-01", tz = "UTC"),
                  length.out = 1000L, by = "2 mins")
  paths <- sim_path_walk(.map = map,
                         .timeline = timeline,
                         .state = "StateXYZD",
                         .xinit = NULL, .n_path = 1L,
                         .model_move = move_xyzd(),
                         .plot = TRUE,
                         .one_page = TRUE)

  #### Simulate observations arising from the simulated path
  # Register a custom `ModelObs` sub-type in Julia
  # * We imagine a pelagic animal in which the depth at each time step
  # * ... is normally distributed around the previous depth.
  # * We write a `ModelObs` sub-type in `Julia` that contains the parameters
  # * ... for this model (i.e., the sigma parameter of the normal distribution).
  julia_command(
    '
    struct ModelObsDepthNormal <: Patter.ModelObs
      sensor_id::Int64
      depth_sigma::Float64
    end
  '
  )
  # Define a `Patter.simulate_obs()` method
  # * We need to specify a function that simulates depths for `ModelObsDepthNormal`
  # * We simulate depths around the previous depth (`state.z`), truncated between
  # * ... the depth of the seabed (`state.map_value`) and the surface.
  julia_command(
    '
  function Patter.simulate_obs(state::StateXYZD, model::ModelObsDepthNormal, t::Int64)
    dbn   = truncated(Normal(state.z, model.depth_sigma), 0, state.map_value)
    rand(dbn)
  end
  '
  )
  # Simulate observations
  obs <- sim_observations(.timeline = timeline,
                          .model_obs = "ModelObsDepthNormal",
                          .model_obs_pars = list(data.table(sensor_id = 1L, depth_sigma = 5)))
  obs <- obs$ModelObsDepthNormal[[1]]
  # Plot simulated depth trajectory
  # * Blue: simulated time series
  # * Grey: seabed depth for simulated time series
  ylim <- range(c(obs$obs, paths$map_value) * -1)
  plot(obs$timestamp, obs$obs * -1, ylim = ylim, col = "royalblue", type = "l")
  lines(paths$timestamp, paths$map_value * -1, col = "grey")

  #### Run the forward filter
  # (optional) Define initial states, by:
  # A) Starting the filter in the correct location by masking `.map`
  # B) Specifying a `map_init()` method based on the observation model
  # C) Specifying a complete data.table of initial state(s)
  origin       <- terra::setValues(map, NA)
  cell         <- terra::cellFromXY(map, cbind(paths$x[1], paths$y[1]))
  origin[cell] <- paths$map_value[1]
  # Define a `Patter.logpdf_obs()` method
  # * This is used to evaluate the log probability of a depth observation
  julia_command(
    '
  function Patter.logpdf_obs(state::State, model::ModelObsDepthNormal, t::Int64, obs::Float64)
    dbn   = truncated(Normal(state.map_value, model.depth_sigma),
                      0.0, state.map_value)
    logpdf(dbn, obs)
  end
  '
  )
  # Run the filter
  fwd <- pf_filter(.map = origin,
                   .timeline = timeline,
                   .state = "StateXYZD",
                   .yobs = list(obs),
                   .model_obs = "ModelObsDepthNormal",
                   .model_move = move_xyzd(),
                   .n_particle = 1000L)
  # Visualise reconstructed time series
  # * Black: particle depths
  # * Blue: simulated time series
  # * Grey: seabed depth for simulated time series
  ylim <- range(c(fwd$states$z, obs$obs, paths$map_value) * -1)
  plot(fwd$states$timestamp, fwd$states$z * -1, ylim = ylim, pch = ".")
  lines(obs$timestamp, obs$obs * -1 , col = "royalblue")
  lines(paths$timestamp, paths$map_value * -1, col = "grey")

}
#> `patter::julia_connect()` called @ 2024-08-01 11:40:06... 
#> ... Running `Julia` setup via `JuliaCall::julia_setup()`... 
#> ... Validating Julia installation... 
#> ... Setting up Julia project... 
#> ... Handling dependencies... 
#> ... `Julia` set up with 8 thread(s). 
#> `patter::julia_connect()` call ended @ 2024-08-01 11:40:14 (duration: ~8 sec(s)). 

#> `patter::pf_filter()` called @ 2024-08-01 11:40:16... 
#> ... 11:40:16: Checking user inputs... 
#> ... 11:40:16: Setting initial states... 
#> Warning: `terra::spatSample()` failed: implementing exhaustive sampling...
#> Warning: [spatSample] fewer samples than requested are available

#> ... 11:40:16: Setting observations... 
#> ... 11:40:16: Running filter... 
#> ... 11:40:18: Collating outputs... 
#> `patter::pf_filter()` call ended @ 2024-08-01 11:40:19 (duration: ~3 sec(s)).