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.SMCType
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.PGType
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_multinomialFunction
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_residualFunction
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_stratifiedFunction
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_systematicFunction
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.ResampleWithESSThresholdType
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.TracedRNGType
TracedRNG{R,N,T}

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

AdvancedPS.splitFunction
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

Internals

Particle Sweep

AdvancedPS.ParticleContainerType

Data 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.