|
| 1 | +# SSMProblems |
| 2 | + |
| 3 | +## Installation |
| 4 | + |
| 5 | +In the `julia` REPL: |
| 6 | + |
| 7 | +```julia |
| 8 | +] add SSMProblems |
| 9 | +``` |
| 10 | + |
| 11 | +## Documentation |
| 12 | + |
| 13 | +`SSMProblems` defines a generic interface for _state space models_ (SSMs). Its |
| 14 | +main objective is to provide a consistent interface for filtering and smoothing |
| 15 | +algorithms to interact with. |
| 16 | + |
| 17 | +Consider a standard (Markovian) state-space model from[^Murray]: |
| 18 | + |
| 19 | + |
| 20 | +[^Murray]: |
| 21 | + > Murray, Lawrence & Lee, Anthony & Jacob, Pierre. (2013). Rethinking resampling in the particle filter on graphics processing units. |
| 22 | +
|
| 23 | +The following three distributions fully specify the model: |
| 24 | + |
| 25 | +- The __initialisation__ distribution, ``f_0``, for the initial latent state ``X_0`` |
| 26 | +- The __transition__ distribution, ``f``, for the latent state ``X_t`` given the previous ``X_{t-1}`` |
| 27 | +- The __observation__ distribution, ``g``, for an observation ``Y_t`` given the state ``X_t`` |
| 28 | + |
| 29 | +The dynamics of the model are given by, |
| 30 | + |
| 31 | +```math |
| 32 | +\begin{aligned} |
| 33 | +x_0 &\sim f_0(x_0) \\ |
| 34 | +x_t | x_{t-1} &\sim f(x_t | x_{t-1}) \\ |
| 35 | +y_t | x_t &\sim g(y_t | x_{t}) |
| 36 | +\end{aligned} |
| 37 | +``` |
| 38 | + |
| 39 | +and the joint law is, |
| 40 | + |
| 41 | +```math |
| 42 | +p(x_{0:T}, y_{0:T}) = f_0(x_0) \prod_t g(y_t | x_t) f(x_t | x_{t-1}). |
| 43 | +``` |
| 44 | + |
| 45 | +We can consider a state space model as being made up of two components: |
| 46 | + |
| 47 | +- A latent Markov chain describing the evolution of the latent state |
| 48 | +- An observation process describing the relationship between the latent states and the observations |
| 49 | + |
| 50 | +Through this lens, we see that the distributions ``f_0``, ``f`` fully describe the latent Markov chain, whereas ``g`` describes the observation process. |
| 51 | + |
| 52 | +A user of `SSMProblems` may define these three distributions directly. |
| 53 | +Alternatively, they can define a subset of methods for sampling and evaluating |
| 54 | +log-densities of the distributions, depending on the requirements of the |
| 55 | +filtering/smoothing algorithms they intend to use. |
| 56 | + |
| 57 | +Using the first approach, we can define a simple linear state space model as follows: |
| 58 | + |
| 59 | +```julia |
| 60 | +using Distributions |
| 61 | +using SSMProblems |
| 62 | + |
| 63 | +struct SimpleLatentDynamics <: LatentDynamics end |
| 64 | + |
| 65 | +function distribution(rng::AbstractRNG, dyn::SimpleLatentDynamics; kwargs...) |
| 66 | + return Normal(0.0, 1.0) |
| 67 | +end |
| 68 | + |
| 69 | +function distribution(rng::AbstractRNG, dyn::SimpleLatentDynamics, step::Int, state::Float64; kwargs...) |
| 70 | + return Normal(state, 0.1) |
| 71 | +end |
| 72 | + |
| 73 | +struct SimpleObservationProcess <: ObservationProcess end |
| 74 | + |
| 75 | +function distribution( |
| 76 | + obs::SimpleObservationPRocess, step::Int, state::Float64, observation::Float64; kwargs... |
| 77 | +) |
| 78 | + return Normal(state, 0.5) |
| 79 | +end |
| 80 | + |
| 81 | +# Construct an SSM from the components |
| 82 | +dyn = SimpleLatentDynamics() |
| 83 | +obs = SimpleObservationProcess() |
| 84 | +model = StateSpaceModel(dyn, obs) |
| 85 | +``` |
| 86 | + |
| 87 | +There are a few things to note here: |
| 88 | + |
| 89 | +- Two methods must be defined for the `LatentDynamics`, one containing |
| 90 | + `step`/`state` arguments and used for transitioning, and one without these, |
| 91 | + used for initialisation. |
| 92 | +- Every function should accept keyword arguments. This is key feature of |
| 93 | + `SSMProblems` that allows it to flexibly represent more exotic models without |
| 94 | + any performance penalty. You can read more about it [here](kwargs.md). |
| 95 | +- If your latent dynamics and observation process cannot be represented as a |
| 96 | + `Distribution` object, you may implement specific methods for sampling and |
| 97 | + log-density evaluation as documented below. |
| 98 | + |
| 99 | +These distribution definitions are used to define `simulate` and `logdensity` |
| 100 | +methods for the latent dynamics and observation process. Package users can then interact with the state space model through these functions. |
| 101 | + |
| 102 | +For example, a bootstrap filter targeting the filtering distribution ``p(x_t | y_{0:t})`` using `N` particles would roughly follow: |
| 103 | + |
| 104 | +```julia |
| 105 | +dyn, obs = model.dyn, model.obs |
| 106 | + |
| 107 | +for (i, observation) in enumerate(observations) |
| 108 | + idx = resample(rng, log_weights) |
| 109 | + particles = particles[idx] |
| 110 | + for i in 1:N |
| 111 | + particles[i] = simulate(rng, dyn, i, particles[i]) |
| 112 | + log_weights[i] += logdensity(obs, i, particles[i], observation) |
| 113 | + end |
| 114 | +end |
| 115 | +``` |
| 116 | + |
| 117 | +For more thorough examples, see the provided example scripts. |
| 118 | + |
| 119 | +### Interface |
| 120 | +```@autodocs |
| 121 | +Modules = [SSMProblems] |
| 122 | +Order = [:type, :function, :module] |
| 123 | +``` |
![github-actions[bot]](https://avatars.githubusercontent.com/in/15368?v=4&size=40){kind=avatar}
0 commit comments