From 885bcfddbd90b1026f1030b5aafb35778d0d44a5 Mon Sep 17 00:00:00 2001 From: audreyolivier Date: Thu, 2 Feb 2023 16:37:16 -0800 Subject: [PATCH] tempering and parallel tempering documentation --- docs/source/bibliography.bib | 17 +++- docs/source/sampling/mcmc/index.rst | 1 + docs/source/sampling/mcmc/tempering.rst | 56 ++++++++++++ .../tempering_mcmc/ParallelTemperingMCMC.py | 91 ++++++++----------- .../tempering_mcmc/baseclass/TemperingMCMC.py | 55 ++++------- 5 files changed, 129 insertions(+), 91 deletions(-) create mode 100644 docs/source/sampling/mcmc/tempering.rst diff --git a/docs/source/bibliography.bib b/docs/source/bibliography.bib index f46fc199a..e01f8d20e 100644 --- a/docs/source/bibliography.bib +++ b/docs/source/bibliography.bib @@ -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}, -} \ No newline at end of file +} + +@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} +} diff --git a/docs/source/sampling/mcmc/index.rst b/docs/source/sampling/mcmc/index.rst index 7ba4f66ae..a54e0f392 100644 --- a/docs/source/sampling/mcmc/index.rst +++ b/docs/source/sampling/mcmc/index.rst @@ -75,6 +75,7 @@ List of MCMC algorithms DRAM DREAM Stretch + Tempering MCMC Adding New MCMC Algorithms diff --git a/docs/source/sampling/mcmc/tempering.rst b/docs/source/sampling/mcmc/tempering.rst new file mode 100644 index 000000000..57073a51e --- /dev/null +++ b/docs/source/sampling/mcmc/tempering.rst @@ -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: diff --git a/src/UQpy/sampling/mcmc/tempering_mcmc/ParallelTemperingMCMC.py b/src/UQpy/sampling/mcmc/tempering_mcmc/ParallelTemperingMCMC.py index 789edb5eb..e9349e596 100644 --- a/src/UQpy/sampling/mcmc/tempering_mcmc/ParallelTemperingMCMC.py +++ b/src/UQpy/sampling/mcmc/tempering_mcmc/ParallelTemperingMCMC.py @@ -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, @@ -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) @@ -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 = [], [] @@ -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: diff --git a/src/UQpy/sampling/mcmc/tempering_mcmc/baseclass/TemperingMCMC.py b/src/UQpy/sampling/mcmc/tempering_mcmc/baseclass/TemperingMCMC.py index 9fba5aafb..ef08bf4df 100644 --- a/src/UQpy/sampling/mcmc/tempering_mcmc/baseclass/TemperingMCMC.py +++ b/src/UQpy/sampling/mcmc/tempering_mcmc/baseclass/TemperingMCMC.py @@ -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