# API

## Samplers

AdvancedPS introduces a few samplers extending AbstractMCMC. The `sample`

method expects a custom type that subtypes `AbstractMCMC.AbstractModel`

. The available samplers are listed below:

### SMC

`AdvancedPS.SMC`

— Type```
SMC(n[, resampler = ResampleWithESSThreshold()])
SMC(n, [resampler = resample_systematic, ]threshold)
```

Create a sequential Monte Carlo (SMC) sampler with `n`

particles.

If the algorithm for the resampling step is not specified explicitly, systematic resampling is performed if the estimated effective sample size per particle drops below 0.5.

The SMC sampler populates a set of particles in a `AdvancedPS.ParticleContainer`

and performs a `AdvancedPS.sweep!`

which propagates the particles and provides an estimation of the log-evidence

```
sampler = SMC(nparticles)
chains = sample(model, sampler)
```

### Particle Gibbs

`AdvancedPS.PG`

— Type```
PG(n, [resampler = ResampleWithESSThreshold()])
PG(n, [resampler = resample_systematic, ]threshold)
```

Create a Particle Gibbs sampler with `n`

particles.

If the algorithm for the resampling step is not specified explicitly, systematic resampling is performed if the estimated effective sample size per particle drops below 0.5.

The Particle Gibbs introduced in ^{[2]} runs a sequence of conditional SMC steps where a pre-selected particle, the reference particle, is replayed and propagated through the SMC step.

```
sampler = PG(nparticles)
chains = sample(model, sampler, nchains)
```

For more detailed examples please refer to the Examples page.

## Resampling

AdvancedPS implements adaptive resampling for both `AdvancedPS.PG`

and `AdvancedPS.SMC`

. The following resampling schemes are implemented:

`AdvancedPS.resample_multinomial`

— Function`resample_multinomial(rng, weights, n)`

Return a vector of `n`

samples `x₁`

, ..., `xₙ`

from the numbers 1, ..., `length(weights)`

, generated by multinomial resampling.

The new indices are sampled from the multinomial distribution with probabilities equal to `weights`

`AdvancedPS.resample_residual`

— Function`resample_residual(rng, weights, n)`

Return a vector of `n`

samples `x₁`

, ..., `xₙ`

from the numbers 1, ..., `length(weights)`

, generated by residual resampling.

In residual resampling we start by duplicating all the particles whose weight is bigger than `1/n`

. We copy each of these particles $N_i$ times where

$N_i = \left\lfloor n w_i \right\rfloor$

We then duplicate the $R_t = n - \sum_i N_i$ missing particles using multinomial resampling with the residual weights given by:

$\tilde{w} = w_i - \frac{N_i}{N}$

`AdvancedPS.resample_stratified`

— Function`resample_stratified(rng, weights, n)`

Return a vector of `n`

samples `x₁`

, ..., `xₙ`

from the numbers 1, ..., `length(weights)`

, generated by stratified resampling.

In stratified resampling `n`

ordered random numbers `u₁`

, ..., `uₙ`

are generated, where

$uₖ \sim U[(k - 1) / n, k / n)$.

Based on these numbers the samples `x₁`

, ..., `xₙ`

are selected according to the multinomial distribution defined by the normalized `weights`

, i.e., `xᵢ = j`

if and only if

$uᵢ \in \left[\sum_{s=1}^{j-1} weights_{s}, \sum_{s=1}^{j} weights_{s}\right)$.

`AdvancedPS.resample_systematic`

— Function`resample_systematic(rng, weights, n)`

Return a vector of `n`

samples `x₁`

, ..., `xₙ`

from the numbers 1, ..., `length(weights)`

, generated by systematic resampling.

In systematic resampling a random number $u \sim U[0, 1)$ is used to generate `n`

ordered numbers `u₁`

, ..., `uₙ`

where

$uₖ = (u + k − 1) / n$.

Based on these numbers the samples `x₁`

, ..., `xₙ`

are selected according to the multinomial distribution defined by the normalized `weights`

, i.e., `xᵢ = j`

if and only if

$uᵢ \in \left[\sum_{s=1}^{j-1} weights_{s}, \sum_{s=1}^{j} weights_{s}\right)$

Each of these schemes is wrapped in a `AdvancedPS.ResampleWithESSThreshold`

struct to trigger a resampling step whenever the ESS is below a certain threshold.

`AdvancedPS.ResampleWithESSThreshold`

— Type`ResampleWithESSThreshold{R,T<:Real}`

Perform resampling using `R`

if the effective sample size is below `T`

. By default we use `resample_systematic`

with a threshold of 0.5

## RNG

AdvancedPS replays the individual trajectories instead of storing the intermediate values. This way we can build efficient samplers. However in order to replay the trajectories we need to reproduce most of the random numbers generated during the execution of the program while also generating diverging traces after each resampling step. To solve these two issues AdvancedPS uses counter-based RNG introduced in ^{[1]} and widely used in large parallel systems see StochasticDifferentialEquations or JAX for other implementations.

Under the hood AdvancedPS is using Random123 for the generators. Using counter-based RNG allows us to split generators thus creating new independent random streams. These generators are also wrapped in a `AdvancedPS.TracedRNG`

type. The `TracedRNG`

keeps track of the keys generated at every `split`

and can be reset to replay random streams.

`AdvancedPS.TracedRNG`

— Type`TracedRNG{R,N,T}`

Wrapped random number generator from Random123 to keep track of random streams during model evaluation

`AdvancedPS.split`

— Function`split(key::Integer, n::Integer=1)`

Split `key`

into `n`

new keys

`AdvancedPS.load_state!`

— Function`load_state!(r::TracedRNG)`

Load state from current model iteration. Random streams are now replayed

`AdvancedPS.save_state!`

— Function`save_state!(r::TracedRNG)`

Add current key of the inner rng in `r`

to `keys`

.

## Internals

### Particle Sweep

`AdvancedPS.ParticleContainer`

— TypeData structure for particle filters

`effectiveSampleSize(pc :: ParticleContainer)`

: Return the effective sample size of the particles in`pc`

`AdvancedPS.sweep!`

— Function`sweep!(rng, pc::ParticleContainer, resampler)`

Perform a particle sweep and return an unbiased estimate of the log evidence.

The resampling steps use the given `resampler`

.

**Reference**

Del Moral, P., Doucet, A., & Jasra, A. (2006). Sequential monte carlo samplers. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 68(3), 411-436.

- 1John K. Salmon, Mark A. Moraes, Ron O. Dror, and David E. Shaw. 2011. Parallel random numbers: as easy as 1, 2, 3. In Proceedings of 2011 International Conference for High Performance Computing, Networking, Storage and Analysis (SC '11). Association for Computing Machinery, New York, NY, USA, Article 16, 1–12. DOI:https://doi.org/10.1145/2063384.2063405
- 2Andrieu, Christophe, Arnaud Doucet, and Roman Holenstein. "Particle Markov chain Monte Carlo methods." Journal of the Royal Statistical Society: Series B (Statistical Methodology) 72, no. 3 (2010): 269-342.