dust provides support for simulating lots of particles at once. Sometimes that will be with a single parameter set and many traces to capture stochastic variability, but other times different particles will need different parameter sets too.
To allow this, dust supports creating a model and indicating that the
parameters represents some group of parameter sets. Consider our SIR
example from vignette("dust"):
sir <- dust::dust_example("sir")For one parameter set with 8 particles we might do:
pars <- list(beta = 0.1)
mod <- sir$new(pars, 0, 8)The “shape” of this model is most easily seen in the state:
mod$state()
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] 1000 1000 1000 1000 1000 1000 1000 1000
#> [2,] 10 10 10 10 10 10 10 10
#> [3,] 0 0 0 0 0 0 0 0
#> [4,] 0 0 0 0 0 0 0 0
#> [5,] 0 0 0 0 0 0 0 0So here we have 5 state variables and 8 particles. The number of state variables is a property of the model, while the number of particles is a property of how you have configured the model, and we will call that the “shape”
mod$shape()
#> [1] 8There are also two helper methods which report the number of particles and number of particles per parameter set, here these will both report 8:
mod$n_particles()
#> [1] 8
mod$n_particles_each()
#> [1] 8So far so easy.
Suppose we have 2 regions (or other grouping) that we want to simulate at once and these have different parameter sets. We create an unnamed list containing the parameter sets
and then when creating the model we indicate that we want 4 particles each for our 2 parameter sets. The total number of particles will be 4 * 2 = 8
mod <- sir$new(pars, 0, 4, pars_multi = TRUE)The state has a different shape now:
mod$state()
#> , , 1
#>
#> [,1] [,2] [,3] [,4]
#> [1,] 1000 1000 1000 1000
#> [2,] 10 10 10 10
#> [3,] 0 0 0 0
#> [4,] 0 0 0 0
#> [5,] 0 0 0 0
#>
#> , , 2
#>
#> [,1] [,2] [,3] [,4]
#> [1,] 1000 1000 1000 1000
#> [2,] 12 12 12 12
#> [3,] 0 0 0 0
#> [4,] 0 0 0 0
#> [5,] 0 0 0 0It’s a 3-dimensional array with dimensions 5, 4, 2 corresponding to
the number of state elements, number of replicates per parameter set and
the number of parameter sets. Note that the second parameter set’s state
is different because the I0 parameter specifies the number
of infected individuals in the initial conditions.
The “shape” of the data here is c(4, 2)
mod$shape()
#> [1] 4 2We still have 8 particles, but 4 per parameter set:
mod$n_particles()
#> [1] 8
mod$n_particles_each()
#> [1] 4This sort of shape is useful where we want to fit multiple particle filters to the data, or simulate replicates blocked across a common set of parameters; that is we’re interested in exploring the stochastic variation that occurs given a particular parameter set.
But there are other way of shaping the state that are useful.
Suppose that you have many parameter sets, each representing (say) an estimate based on some inference process
For this we might write:
mod <- sir$new(pars, 0, 1, pars_multi = TRUE)to indicate 1 particle per parameter. However, that’s not ideal:
dim(mod$state())
#> [1] 5 1 8
mod$shape()
#> [1] 1 8By doing this we end up with a weird dimension in the middle of our state. This would make sense if we were going to also replicate within the parameter set but does not make any sense this way.
To indicate that we don’t want this, pass NULL in for
the number of particles:
mod <- sir$new(pars, 0, NULL, pars_multi = TRUE)
mod$state()
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] 1000 1000 1000 1000 1000 1000 1000 1000
#> [2,] 10 10 10 10 10 10 10 10
#> [3,] 0 0 0 0 0 0 0 0
#> [4,] 0 0 0 0 0 0 0 0
#> [5,] 0 0 0 0 0 0 0 0
mod$shape()
#> [1] 8(why NULL? because if you tried c(5, 1, 8)
and c(5, NULL, 8) you would see the same values as the two
state vectors above).
Further options are useful! Suppose in our previous case, the 8 parameter sets were themselves structured in some block, representing two groups of four parameter sets. To indicate this we create a list with dimensions:
dim(pars) <- c(4, 2)
pars
#> [,1] [,2]
#> [1,] list,1 list,1
#> [2,] list,1 list,1
#> [3,] list,1 list,1
#> [4,] list,1 list,1Then, when creating the object these dimensions will be applied to the data
mod <- sir$new(pars, 0, NULL, pars_multi = TRUE)
mod$state()
#> , , 1
#>
#> [,1] [,2] [,3] [,4]
#> [1,] 1000 1000 1000 1000
#> [2,] 10 10 10 10
#> [3,] 0 0 0 0
#> [4,] 0 0 0 0
#> [5,] 0 0 0 0
#>
#> , , 2
#>
#> [,1] [,2] [,3] [,4]
#> [1,] 1000 1000 1000 1000
#> [2,] 10 10 10 10
#> [3,] 0 0 0 0
#> [4,] 0 0 0 0
#> [5,] 0 0 0 0
mod$shape()
#> [1] 4 2These options have implications for a number of methods, and for the basic running of the simulation. For example, when we run our model with structured parameters:
res <- mod$simulate(0:10)
dim(res)
#> [1] 5 4 2 11Here the simulation output has dimensions:
mod$shape(), as the first
dimension of the parametersmod$shape(), as the second
dimension of the parametersConversely if we’d also replicated parameters with this set of parameters we might see
mod <- sir$new(pars, 0, 3, pars_multi = TRUE)
mod$shape()
#> [1] 3 4 2
dim(mod$simulate(0:10))
#> [1] 5 3 4 2 11Again here note that mod$shape() is sandwiched between
the number of state variables and number of timesteps.
There are a few helpful things on single parameter objects that become ambiguous or awkward with multi-parameter objects. Starting with a single parameter object:
pars <- list(beta = 0.1)
mod <- sir$new(pars, 0, 8)We might set all particles to have the same initial state:
mod$update_state(state = c(1000, 10, 0, 0, 0))
mod$state()
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] 1000 1000 1000 1000 1000 1000 1000 1000
#> [2,] 10 10 10 10 10 10 10 10
#> [3,] 0 0 0 0 0 0 0 0
#> [4,] 0 0 0 0 0 0 0 0
#> [5,] 0 0 0 0 0 0 0 0or we might set each particle to a different state, perhaps representing some stochasticity in starting values
n <- rpois(8, 10)
mod$update_state(state = rbind(1010 - n, n, 0, 0, 0))
mod$state()
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] 1001 993 998 1003 1002 1008 1002 1000
#> [2,] 9 17 12 7 8 2 8 10
#> [3,] 0 0 0 0 0 0 0 0
#> [4,] 0 0 0 0 0 0 0 0
#> [5,] 0 0 0 0 0 0 0 0The difference here is the dimension of the new state. Recall that the state will always be in the form
n_state, shape[1], ..., shape[n]If the dimension lacks the first element of shape, and if there is more than one particle per parameter set, then we replicate state across all particles that share that parameter set.