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
— TypeSMC(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
— TypePG(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
— Functionresample_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
— Functionresample_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
— Functionresample_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
— Functionresample_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
— TypeResampleWithESSThreshold{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
— TypeTracedRNG{R,N,T}
Wrapped random number generator from Random123 to keep track of random streams during model evaluation
AdvancedPS.split
— Functionsplit(key::Integer, n::Integer=1)
Split key
into n
new keys
AdvancedPS.load_state!
— Functionload_state!(r::TracedRNG)
Load state from current model iteration. Random streams are now replayed
AdvancedPS.save_state!
— Functionsave_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 inpc
AdvancedPS.sweep!
— Functionsweep!(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.