Sampling¶

This section explores the various options that are available to the user when it comes to the initialization of the sampler and running the sampler.

Sampler Initialization¶

We assume that the user has already defined their prior probability prior and log-likelihood function likelihood, and if the latter requires any additional arguments (e.g., vectorize, likelihood_args, likelihood_kwargs) to provide to the sampler the user has already done it. Please see the relevant sections in the documentation if you have not done so already.

Number of Particles¶

By far, the most important parameters that control the performance of pocoMC are the number of effective and active particles. These are controlled respectively by the parameters n_effective and n_active.

Number of Effective Particles (n_effective): Monte Carlo samplers such as pocoMC represent probability distributions as ensembles of particles. The greater the number of particles, the better the more accurate the representation of the target distribution. Generally, simpler distributions (i.e., close to Gaussian) require low number of particles, whereas complex distributions (i.e., skewed or multimodal) require large number of particles. Unlike Nested Sampling methods, Preconditioned Monte Carlo, and thus pocoMC, utilize unequally weighted particles with some particles contributing more and some less to the aforementioned representation. As a result, we talk about effective particles. We recommend values around 500 for most problems and increasing that to 2000-4000 particularly challenging distributions. If you are not sure about what value to use, do not worry. This can be determined automatically by pocoMC (see the Dynamic Particle Allocation section below). The default value is n_effective=512.
Number of Active Particles (n_active): pocoMC propages the effective particles from the prior to the posterior through a series of steps or iterations. However, in each iteration, not all effective particles are updated. Instead, only a subset of those, called active particles are updated. Generally, the number of active particles n_active should be no more than half the number of effective particles n_effective for best performance. The default value is n_active=256.

Here is an example of how to manually specify the number of effective and active particles. If only one of them is specified and the other is None, then the other is set automatically.

sampler = pc.Sampler(prior, likelihood, n_effective=1024, n_active=512)

Dynamic Particle Allocation¶

For most applications, the default number of effective and active particles are more than sufficient to guarantee fast and accurate sampling of the posterior distribution. However, there are cases where the difficulty of the target distribution varies during the run. In this case it is beneficial to dynamically adjust the number of effective particles n_effective automatically (i.e., increase it when the problem becomes more difficult and decrease it when it becomes easier). This dynamic allocation of effective particles is achieved by setting dynamic=True in the initialization of the sampler. By default, dynamic allocation is turned on.

sampler = pc.Sampler(prior, likelihood, dynamic=True)

Markov chain Monte Carlo¶

pocoMC relies on Markov chain Monte Carlo (MCMC) methods to diversify the active particles in each iteration. This is the most computationally expensive part of the algorithm and great care has been taken so that it is performed as efficiently as possible.

pocoMC supports two different Markov kernels that can be specified using the sample argument:

t-preconditioned Crank-Nicolson (sample='tpcn'): The t-preconditioned Crank-Nicolson (tpCN) MCMC algorithm is an advanced technique designed to improve sampling efficiency in high-dimensional and complex distributions. Unlike traditional methods (e.g, Random-walk Metropolis) where Gaussian perturbations are added directly to the current state, tpCN shifts these perturbations toward higher probability regions, enhancing the algorithm’s scalability to high dimensions. This shift is achieved by combining the current state with a preconditioned perturbation, balancing exploration and stability. The preconditioning step adapts to the structure of the target distribution, ensuring consistent acceptance rates and efficient exploration. Consequently, tpCN provides a robust and effective approach for sampling in high-dimensional spaces, outperforming simpler MCMC methods in such challenging settings.
Random-walk Metropolis (sample='rwm'): The is an MCMC method designed to sample from complex probability distributions. It generates new samples by adding a Gaussian-distributed perturbation to the current state and then accepts or rejects the proposed state based on the Metropolis acceptance criterion, which ensures that the samples converge to the target distribution. Despite its simplicity and broad applicability, the Random-walk Metropolis algorithm scales poorly to high-dimensional spaces. As the dimensionality increases, the probability of accepting proposed moves decreases, leading to inefficient exploration of the target distribution and slow convergence. This makes it less suitable for high-dimensional problems where more sophisticated methods (e.g., t-preconditioned Crank-Nicolson, NUTS, etc.) are often required.

By default, sample='tpcn' and the t-preconditioned Crank-Nicolson Markov kernel is utilized. The Random-walk Metropolis kernel is only intended to be used for testing.

The purpose of the Markov kernel is to diversify the particles and make the equilibrate in each target distribution. This means that in each iteration we need to ensure that MCMC has ran long enough. To do this, pocoMC monitors the mean (unnormalized) posterior log-probability until it stops increasing for at least N steps, where N is given by:

\[ N = n_{\text{steps}}\times\left(\frac{2.38/\sqrt{n_{\text{dim}}}}{\sigma}\right)^{2} \]

where n_dim is the dimensionality of the problem (i.e., number of parameteres) and \(\sigma\) is the proposal scale of tpCN or RWM that is determined adaptively by pocoMC and n_steps is a hyperparameter. By default, it is equal to the dimensionality of the problem, that is n_steps=n_dim. However, for highly challenging problems, the user can provide a different value. For instance:

sampler = pc.Sampler(prior, likelihood, n_steps=2*n_dim)

Preconditioning¶

By default, pocoMC utilizes a normalizing flow to precondition (i.e., simplify) the target distribution. However, for sufficiently nice distributions that do not deviate significantly from a Gaussian distribution, this procedure may be unnecessary. In these cases, when the posterior is close to a Gaussian distribution and the computational cost of the likelihood is low, one can turn off the preconditioning and avoid training and evaluating the normalizing flow. This can be achieve simply by doing:

sampler = pc.Sampler(prior, likelihood, precondition=False)

Running¶

Once the sampler has been configured and initialized, running it in order to sample from the posterior distribution is quite simple. All the user has to do is call the run() method of the Sampler object. Optionally, there are three parameters that can be adjusted:

n_total : This parameter specifies the effective sample size (or unique sample size if metric='uss') that is required before the sampling terminates. The default value is n_total=4096. This is usually more than enough to construct nice diagrams showing the 1D and 2D marginal posteriors of various parameters. However, for publication purposes one may want to increase n_total to about ~ 10_000. Similarly, for preliminary runs one can also decrese this (e.g., n_total=1024).
n_evidence : This parameter specified the number of samples that will be generated from the normalizing flow at the end of the run to compute the importance sampling estimate of the log model evidence \(\log\mathcal{Z}\). Higher values will generally lead to more accurate estimates. If you do not require an estimate of the model evidence, you can set n_evidence=0.
progress : This is a simple boolean parameter controlling whether or not a progress bar appears during running. By default progress=True.

sampler.run(
    n_total=4096,
    n_evidence=4096,
    progress=True,
)

Continue Running after Completion¶

It is possible to continue running the sampler after sampling has been completed in order to add more samples. This can be useful if the user requires more samples to be able to approximate posterior or estimate the evidence more accurately. This can be achieved easily by calling the run() method again with higher n_total and/or n_evidence values. For instance:

sampler.run(
    n_total=16384,
    n_evidence=16384,
)