Skip to content

Commit

Permalink
tempering and parallel tempering documentation
Browse files Browse the repository at this point in the history
  • Loading branch information
AudOlivier committed Feb 3, 2023
1 parent 5e9dbed commit 885bcfd
Show file tree
Hide file tree
Showing 5 changed files with 129 additions and 91 deletions.
17 changes: 16 additions & 1 deletion docs/source/bibliography.bib
Original file line number Diff line number Diff line change
Expand Up @@ -827,4 +827,19 @@ @article{saltelli_2002
url = {https://www.sciencedirect.com/science/article/pii/S0010465502002801},
author = {Andrea Saltelli},
keywords = {Sensitivity analysis, Sensitivity measures, Sensitivity indices, Importance measures},
}
}

@article{PTMCMC1,
title = {Parallel Tempering: Theory, Applications, and New Perspectives},
author = {David J. Earl, Michael W. Deem},
year = {2005},
doi = {https://doi.org/10.48550/arXiv.physics/0508111}
}

@inproceedings{PTMCMC2,
title = {Using Thermodynamic Integration to Calculate the Posterior Probability in Bayesian Model Selection Problems},
booktitle = {AIP Conference Proceedings 707, 59},
year = {2004},
doi = {https://doi.org/10.1063/1.1751356},
author = {Paul M. Goggans and Ying Chi}
}
1 change: 1 addition & 0 deletions docs/source/sampling/mcmc/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -75,6 +75,7 @@ List of MCMC algorithms
DRAM <dram>
DREAM <dream>
Stretch <stretch>
Tempering MCMC <tempering>


Adding New MCMC Algorithms
Expand Down
56 changes: 56 additions & 0 deletions docs/source/sampling/mcmc/tempering.rst
Original file line number Diff line number Diff line change
@@ -0,0 +1,56 @@
Tempering MCMC
~~~~~~~~~~~~~~

Tempering MCMC algorithms aim at sampling from a target distribution :math:`p_1(x)` of the form
:math:`p_1(x)=\frac{q_1(x)p_0(x)}{Z_1}` where the factor :math:`q_1(x)` and reference distribution
:math:`p_0(x)` can be evaluated. Additionally, these algorithms return an estimate of the normalization
constant :math:`Z_1=\int{q_{1}(x) p_{0}(x)dx}`.

The algorithms sample from a sequence of intermediate densities
:math:`p_{\beta}(x) \propto q_{\beta}(x) p_{0}(x)` for values of the parameter :math:`\beta` between 0 and 1
(:math:`\beta=\frac{1}{T}` where :math:`T` is sometimes called the temperature, :math:`q_{\beta}(x)` is referred to as the intermediate factor associated with tempering parameter :math:`\beta`).
Setting :math:`\beta = 1` equates sampling from the target, while :math:`\beta \rightarrow 0` samples from the
reference distribution.

Parallel tempering samples from all distributions simultaneously, and the tempering parameters :math:`0 < \beta_1 < \beta_2 < \cdots < \beta_{N} \leq 1` must be
chosen in advance by the user. Sequential tempering on the other hand samples from the various distributions sequentially, starting
from the reference distribution, and the tempering parameters are selected adaptively by the algorithm.

The :class:`.TemperingMCMC` base class defines inputs that are common to parallel and sequential tempering:

.. autoclass:: UQpy.sampling.mcmc.tempering_mcmc.TemperingMCMC
:members:
:exclude-members: run, evaluate_normalization_constant

Parallel Tempering
^^^^^^^^^^^^^^^^^^^^

This algorithm (see e.g. :cite:`PTMCMC1` for theory about parallel tempering) runs the chains sampling from the various tempered distributions simultaneously. Periodically during the
run, the different temperatures swap members of their ensemble in a way that preserves detailed balance. The chains
closer to the reference chain (hot chains) can sample from regions that have low probability under the target and
thus allow a better exploration of the parameter space, while the cold chains can better explore regions of high
likelihood.

In parallel tempering, the normalizing constant :math:`Z_1` is evaluated via thermodynamic integration (:cite:`PTMCMC2`). Defining
the potential function :math:`U_{\beta}(x)=\frac{\partial \log{q_{\beta}(x)}}{\partial \beta}`, then

:math:`\ln{Z_1} = \log{Z_0} + \int_{0}^{1} E_{x \sim p_{\beta}} \left[ U_{\beta}(x) \right] d\beta`

where the expectations are approximated via MC sampling using saved samples from the intermediate distributions. A trapezoidal rule is used for integration.

The :class:`.ParallelTemperingMCMC` class is imported using the following command:

>>> from UQpy.sampling.mcmc.tempering_mcmc.ParallelTemperingMCMC import ParallelTemperingMCMC

.. autoclass:: UQpy.sampling.mcmc.tempering_mcmc.ParallelTemperingMCMC
:members:

Sequential Tempering
^^^^^^^^^^^^^^^^^^^^

The :class:`.SequentialTemperingMCMC` class is imported using the following command:

>>> from UQpy.sampling.mcmc.tempering_mcmc.SequentialTemperingMCMC import SequentialTemperingMCMC

.. autoclass:: UQpy.sampling.mcmc.tempering_mcmc.SequentialTemperingMCMC
:members:
91 changes: 36 additions & 55 deletions src/UQpy/sampling/mcmc/tempering_mcmc/ParallelTemperingMCMC.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,42 +7,6 @@


class ParallelTemperingMCMC(TemperingMCMC):
"""
Parallel-Tempering MCMC
This algorithm runs the chains sampling from various tempered distributions in parallel. Periodically during the
run, the different temperatures swap members of their ensemble in a way that preserves detailed balance.The chains
closer to the reference chain (hot chains) can sample from regions that have low probability under the target and
thus allow a better exploration of the parameter space, while the cold chains can better explore the regions of high
likelihood.
In parallel tempering, the normalizing constant :math:`Z_1` is evaluated via thermodynamic integration. Define
the potential function :math:`U_{\beta}(x)=\frac{\partial \log{q_{\beta}(x)}}{\partial \beta}`, then
:math:`\log{Z_{1}} = \log{Z_{0}} + \int_{0}^{1} E_{x~p_{beta}} \left[ U_{\beta}(x) \right] d\beta`
where the expectations are approximated via MC sampling using samples from the intermediate distributions (see
:method:`evaluate_normalization_constant`).
**References**
1. Parallel Tempering: Theory, Applications, and New Perspectives, Earl and Deem
2. Adaptive Parallel Tempering MCMC
3. emcee the MCMC Hammer python package
**Inputs:**
Many inputs are similar to MCMC algorithms (`n_samples`, `n_samples_per_chain`, 'random_state')
:param save_log_pdf: boolean, see :class:MCMC documentation. Importantly, this needs to be set to True if one wants
to evaluate the normalization constant via thermodynamic integration.
:param n_iterations_between_sweeps: number of iterations (sampling steps) between sweeps between chains
:param tempering_parameters: list of :math:`\beta` values so that
:math:`0 < \beta_1 < \beta_2 < \cdots < \beta_N \leq 1`. Either `tempering_parameters` or `n_tempering_parameters`
should be provided.
:param n_tempering_parameters: number of tempering levels N, the tempering parameters are selected to follow a
geometric suite :math:`\frac{1}{\sqrt{2}^(N-n)}` for n in 1:N.
:param samplers: :class:`MCMC` object or list of such objects: MCMC samplers used to sample the parallel chains. If
only one object is provided, the same MCMC sampler is used for all chains. Default to running a simple MH algorithm,
where the proposal covariance for a given chain is :math:`\frac{1}{\beta}`.
"""

@beartype
def __init__(self, n_iterations_between_sweeps: PositiveInteger,
Expand All @@ -55,6 +19,27 @@ def __init__(self, n_iterations_between_sweeps: PositiveInteger,
n_tempering_parameters: int = None,
samplers: Union[MCMC, list[MCMC]] = None):

"""
Class for Parallel-Tempering MCMC.
:param save_log_pdf: boolean, see :class:`MCMC` documentation. Importantly, this needs to be set to True if
one wants to evaluate the normalization constant via thermodynamic integration.
:param n_iterations_between_sweeps: number of iterations (sampling steps) between sweeps between chains.
:param tempering_parameters: tempering parameters, as a list of N floats increasing from 0. to 1. Either
`tempering_parameters` or `n_tempering_parameters` should be provided
:param n_tempering_parameters: number of tempering levels N, the tempering parameters are selected to follow
a geometric suite by default
:param samplers: :class:`MCMC` object or list of such objects: MCMC samplers used to sample the parallel
chains. If only one object is provided, the same MCMC sampler is used for all chains. Default to running a
simple MH algorithm, where the proposal covariance for a given chain is inversely proportional to the
tempering parameter.
"""

super().__init__(pdf_intermediate=pdf_intermediate, log_pdf_intermediate=log_pdf_intermediate,
args_pdf_intermediate=args_pdf_intermediate, distribution_reference=None,
save_log_pdf=save_log_pdf, random_state=random_state)
Expand Down Expand Up @@ -121,14 +106,14 @@ def run(self, nsamples: PositiveInteger = None, nsamples_per_chain: PositiveInte
"""
Run the MCMC algorithm.
This function samples from the MCMC chains and appends samples to existing ones (if any). This method leverages
the :method:`run_iterations` method specific to each of the samplers.
This function samples from the MCMC chains and appends samples to existing ones (if any). This method
leverages the `run_iterations` method specific to each of the samplers.
:param nsamples: Number of samples to generate from the target (the same number of samples will be generated
for all intermediate distributions).
:param nsamples_per_chain: Number of samples to generate per chain of the target MCMC sampler. Either `nsamples`
or `nsamples_per_chain` must be provided (not both).
:param nsamples_per_chain: Number of samples per chain to generate from the target. Either
`nsamples` or `nsamples_per_chain` must be provided (not both)
"""
current_state, current_log_pdf = [], []
Expand Down Expand Up @@ -204,23 +189,19 @@ def run(self, nsamples: PositiveInteger = None, nsamples_per_chain: PositiveInte
@beartype
def evaluate_normalization_constant(self, compute_potential, log_Z0: float = None, nsamples_from_p0: int = None):
"""
Evaluate normalization constant :math:`Z_1` as:
:math:`\log{Z_{1}} \approx \log{Z_{\beta_{1}}} + \int_{\beta_{1}}^{\beta_{N}} \frac{1}{M} \sum_m U_{\beta}(x_m)
d\beta`
where :math:`x_m` are samples from the intermediate distributions. The integration is performed via a
trapezoidal rule. The function returns an approximation of :math:`Z_1`, and also saves intermediate results
(value of :math:`log_Z0`, list of tempering parameters used in integration, and values of the associated
expected potentials :math:`E_{x \sim p_{\beta} \left[ U_{\beta}(x) \right]}\frac{1}{M} \sum_m U_{\beta}(x_m)`)
Evaluate normalization constant :math:`Z_1`.
The function returns an approximation of :math:`Z_1`, and saves intermediate results
(value of :math:`\ln{Z_0}`, list of tempering parameters used in integration, and values of the associated
expected potentials.
:param compute_potential: Function that takes three inputs: :code:`x` (sample points where to evaluate the
potential), :code:`log_factor_tempered_values` (values of :math:`\log{q_{\beta}(x)}` at these points), `beta`
(tempering parameter) and evaluates the potential
:math:`U_{\beta}(x)=\frac{\partial \log{q_{\beta}(x)}}{\partial \beta}`.
:param log_Z0: Value of :math:`\log{Z_{0}}` (float), if unknwon, see `nsamples_from_p0`
:param nsamples_from_p0: number of samples :math:`M_0` from the reference distribution :math:`p_0` to be used
to evaluate :math:`\log{Z_{0}} \approx \frac{1}{M_0} \sum{p_{\beta=0}(x)}`. Used only if input `log_Z0` is not
provided.
potential), :code:`log_factor_tempered_values` (values of the log intermediate factors evaluated at points
:code:`x`), :code:`temper_param`(tempering parameter) and evaluates the potential:
:param log_Z0: Value of :math:`\ln{Z_{0}}` (float), if unknwon, see `nsamples_from_p0`.
:param nsamples_from_p0: number of samples from the reference distribution to sample to evaluate :math:`\ln{Z_{0}}`. Used only if input `log_Z0` is not provided.
"""
if not self.save_log_pdf:
Expand Down
55 changes: 20 additions & 35 deletions src/UQpy/sampling/mcmc/tempering_mcmc/baseclass/TemperingMCMC.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,44 +3,29 @@


class TemperingMCMC(ABC):
"""
Parent class to parallel and sequential tempering MCMC algorithms.
These algorithms aim at sampling from a target distribution :math:`p_1(x)` of the form
:math:`p_1(x)=\frac{q_1(x)p_0(x)}{Z_1}` where the intermediate factor :math:`q_1(x)` and reference distribution
:math:`p_0(x)` can be evaluated. Additionally, these algorithms return an estimate of the normalization
constant :math:`Z_1=\int{q_{1}(x) p_{0}(x)dx}`.
The algorithms sample from a sequence of intermediate densities
:math:`p_{\beta}(x) \propto q_{\beta}(x) p_{0}(x)` for values of the parameter :math:`\beta` between 0 and 1
(:math:`\beta=\frac{1}{T}` where :math:`T` is sometimes called the temperature).
Setting :math:`\beta = 1` equates sampling from the target, while :math:`\beta \rightarrow 0` samples from the
reference distribution.
Parallel tempering samples from all distributions simultaneously, and the tempering parameters :math:`\beta` must be
chosen in advance by the user. Sequential tempering samples from the various distributions sequentially, starting
from the reference distribution, and the tempering parameters are selected adaptively by the algorithm.
Inputs that are common to both parallel and sequential tempering algorithms are:
:param pdf_intermediate: callable that computes the intermediate factor :math:`q_{\beta}(x)`. It should take at
least two inputs :code:`x` (ndarray, point(s) at which to evaluate the function), and :code:`temper_param` (float,
tempering parameter :math:`\beta`). Either `pdf_intermediate` or `log_pdf_intermediate` must be provided
(`log_pdf_intermediate` is preferred). Within the code, the `log_pdf_intermediate` is evaluated as:
:code:`log q_{\beta}(x) = log_pdf_intermediate(x, \beta, *args_pdf_intermediate)`
where `args_pdf_intermediate` are additional positional arguments that are provided to the class via its
`args_pdf_intermediate` input.
:param log_pdf_intermediate: see `pdf_intermediate`
:param args_pdf_intermediate: see `pdf_intermediate`
:param distribution_reference: reference pdf :math:`p_0` as a :class:`.Distribution` object
:param save_log_pdf: see same input in :class:`MCMC`
:param random_state
"""

def __init__(self, pdf_intermediate=None, log_pdf_intermediate=None, args_pdf_intermediate=(),
distribution_reference=None, save_log_pdf=True, random_state=None):
"""
Parent class to parallel and sequential tempering MCMC algorithms.
:param pdf_intermediate: callable that computes the intermediate factor. It should take at
least two inputs :code:`x` (ndarray, point(s) at which to evaluate the function), and :code:`temper_param` (float,
tempering parameter). Either `pdf_intermediate` or `log_pdf_intermediate` must be provided
(`log_pdf_intermediate` is preferred). Within the code, the `log_pdf_intermediate` is evaluated as:
:code:`log_pdf_intermediate(x, temper_param, *args_pdf_intermediate)`
where `args_pdf_intermediate` are additional positional arguments that are provided to the class via its
`args_pdf_intermediate` input
:param log_pdf_intermediate: see `pdf_intermediate`
:param args_pdf_intermediate: see `pdf_intermediate`
:param distribution_reference: reference pdf :math:`p_0` as a :class:`.Distribution` object
:param save_log_pdf: see same input in :class:`MCMC`
"""
self.logger = logging.getLogger(__name__)
# Check a few inputs
self.save_log_pdf = save_log_pdf
Expand Down

0 comments on commit 885bcfd

Please sign in to comment.