Rao-Blackwell adaptive MCMC

docs/source/mcmc_samplers/index.rst

@@ -25,6 +25,7 @@ interface, that can be used to sample from an unknown
     metropolis_mcmc
     monomial_gamma_hamiltonian_mcmc
     population_mcmc
+    rao_blackwell_ac_mcmc
     relativistic_mcmc
     slice_doubling_mcmc
     slice_stepout_mcmc

docs/source/mcmc_samplers/rao_blackwell_ac_mcmc.rst

@@ -0,0 +1,7 @@
+******************
+Rao-Blackwell ACMC
+******************
+
+.. currentmodule:: pints
+
+.. autoclass:: RaoBlackwellACMC

examples/README.md

@@ -44,6 +44,7 @@ relevant code.
 - [Haario-Bardenet Adaptive Covariance MCMC](./sampling-adaptive-covariance-haario-bardenet.ipynb)
 - [Metropolis Random Walk MCMC](./sampling-metropolis-mcmc.ipynb)
 - [Population MCMC](./sampling-population-mcmc.ipynb)
+- [Rao-Blackwell Adaptive Covariance MCMC](./sampling-adaptive-covariance-rao-blackwell.ipynb)
 - [Slice Sampling: Stepout MCMC](./sampling-slice-stepout-mcmc.ipynb)
 - [Slice Sampling: Doubling MCMC](./sampling-slice-doubling-mcmc.ipynb)
 - [Slice Sampling: Overrelaxation MCMC](./sampling-slice-overrelaxation-mcmc.ipynb)

pints/__init__.py

@@ -208,6 +208,7 @@ def version(formatted=False):
 from ._mcmc._metropolis import MetropolisRandomWalkMCMC
 from ._mcmc._monomial_gamma_hamiltonian import MonomialGammaHamiltonianMCMC
 from ._mcmc._population import PopulationMCMC
+from ._mcmc._rao_blackwell_ac import RaoBlackwellACMC
 from ._mcmc._relativistic import RelativisticMCMC
 from ._mcmc._slice_stepout import SliceStepoutMCMC
 from ._mcmc._slice_doubling import SliceDoublingMCMC

pints/_mcmc/_rao_blackwell_ac.py

@@ -0,0 +1,110 @@
+#
+# Rao-Blackwell Adaptive MCMC method
+#
+# This file is part of PINTS.
+#  Copyright (c) 2017-2019, University of Oxford.
+#  For licensing information, see the LICENSE file distributed with the PINTS
+#  software package.
+#
+from __future__ import absolute_import, division
+from __future__ import print_function, unicode_literals
+import pints
+import numpy as np
+
+
+class RaoBlackwellACMC(pints.GlobalAdaptiveCovarianceMC):
+    """
+    Rao-Blackwell adaptive MCMC, as described by Algorithm 3 in [1]_.
+    After initialising mu0 and sigma0, in each iteration after initial
+    phase (t), the following steps occur::
+
+        theta* ~ N(theta_t, lambda * sigma0)
+        alpha(theta_t, theta*) = min(1, p(theta*|data) / p(theta_t|data))
+        u ~ uniform(0, 1)
+        if alpha(theta_t, theta*) > u:
+            theta_t+1 = theta*
+        else:
+            theta_t+1 = theta_t
+        mu_t+1 = mu_t + gamma_t+1 * (theta_t+1 - mu_t)
+        sigma_t+1 = sigma_t + gamma_t+1 *
+                        (bar((theta_t+1 - mu_t)(theta_t+1 - mu_t)') - sigma_t)
+
+    where::
+
+        bar(theta_t+1) = alpha(theta_t, theta*) theta* +
+                            (1 - alpha(theta_t, theta*)) theta_t
+
+    Note that we deviate from the paper in two places::
+
+        gamma_t = t^-eta
+        Y_t+1 ~ N(theta_t, lambda * sigma0) rather than
+            Y_t+1 ~ N(theta_t, sigma0)
+
+    Extends :class:`GlobalAdaptiveCovarianceMC`.
+
+    References
+    ----------
+    .. [1] A tutorial on adaptive MCMC
+           Christophe Andrieu and Johannes Thoms, Statistical Computing, 2008,
+           18: 343-373.
+           https://doi.org/10.1007/s11222-008-9110-y
+    """
+    def __init__(self, x0, sigma0=None):
+        super(RaoBlackwellACMC, self).__init__(x0, sigma0)
+
+        # heuristic based on normal approximation
+        self._lambda = (2.38**2) / self._n_parameters
+
+        self._X = None
+        self._Y = None
+
+    def ask(self):
+        """ See :meth:`SingleChainMCMC.ask()`. """
+        super(RaoBlackwellACMC, self).ask()
+
+        # Propose new point
+        if self._proposed is None:
+            self._proposed = np.random.multivariate_normal(
+                self._current, self._lambda * self._sigma)
+
+            # Set as read-only
+            self._proposed.setflags(write=False)
+
+        # Return proposed point
+        return self._proposed
+
+    def n_hyper_parameters(self):
+        """ See :meth:`TunableMethod.n_hyper_parameters()`. """
+        return 1
+
+    def name(self):
+        """ See :meth:`pints.MCMCSampler.name()`. """
+        return 'Rao-Blackwell adaptive covariance MCMC'
+
+    def tell(self, fx):
+        """ See :meth:`pints.AdaptiveCovarianceMCMC.tell()`. """
+        self._Y = np.copy(self._proposed)
+        self._X = np.copy(self._current)
+
+        super(RaoBlackwellACMC, self).tell(fx)
+
+        return self._current
+
+    def _update_sigma(self):
+        """
+        Updates sigma using Rao-Blackwellised formula::
+
+            sigma_t+1 = sigma_t + gamma_t+1 *
+                        (bar((theta_t+1 - mu_t)(theta_t+1 - mu_t)') - sigma_t)
+
+        where::
+
+            bar(X_t+1) = alpha(X_t, Y_t+1) * Y_t+1 +
+                            (1 - alpha(X_t, Y_t+1)) * X_t
+        """
+        acceptance_prob = (
+            np.minimum(1, np.exp(self._log_acceptance_ratio)))
+        X_bar = acceptance_prob * self._Y + (1.0 - acceptance_prob) * self._X
+        dsigm = np.reshape(X_bar - self._mu, (self._n_parameters, 1))
+        self._sigma = ((1 - self._gamma) * self._sigma +
+                       self._gamma * np.dot(dsigm, dsigm.T))

pints/tests/test_mcmc_rao_blackwell_ac.py

@@ -0,0 +1,129 @@
+#!/usr/bin/env python
+#
+# Tests the basic methods of the Rao-Blackwel adaptive covariance MCMC routine.
+#
+# This file is part of PINTS.
+#  Copyright (c) 2017-2019, University of Oxford.
+#  For licensing information, see the LICENSE file distributed with the PINTS
+#  software package.
+#
+import pints
+import pints.toy as toy
+import unittest
+import numpy as np
+
+from shared import StreamCapture
+
+# Consistent unit testing in Python 2 and 3
+try:
+    unittest.TestCase.assertRaisesRegex
+except AttributeError:
+    unittest.TestCase.assertRaisesRegex = unittest.TestCase.assertRaisesRegexp
+
+
+class TestRaoBlackwellACMC(unittest.TestCase):
+    """
+    Tests the basic methods of the Rao-Blackwell AC routine.
+    """
+
+    @classmethod
+    def setUpClass(cls):
+        """ Set up problem for tests. """
+
+        # Create toy model
+        cls.model = toy.LogisticModel()
+        cls.real_parameters = [0.015, 500]
+        cls.times = np.linspace(0, 1000, 1000)
+        cls.values = cls.model.simulate(cls.real_parameters, cls.times)
+
+        # Add noise
+        cls.noise = 10
+        cls.values += np.random.normal(0, cls.noise, cls.values.shape)
+        cls.real_parameters.append(cls.noise)
+        cls.real_parameters = np.array(cls.real_parameters)
+
+        # Create an object with links to the model and time series
+        cls.problem = pints.SingleOutputProblem(
+            cls.model, cls.times, cls.values)
+
+        # Create a uniform prior over both the parameters and the new noise
+        # variable
+        cls.log_prior = pints.UniformLogPrior(
+            [0.01, 400, cls.noise * 0.1],
+            [0.02, 600, cls.noise * 100]
+        )
+
+        # Create a log likelihood
+        cls.log_likelihood = pints.GaussianLogLikelihood(cls.problem)
+
+        # Create an un-normalised log-posterior (log-likelihood + log-prior)
+        cls.log_posterior = pints.LogPosterior(
+            cls.log_likelihood, cls.log_prior)
+
+    def test_method(self):
+
+        # Create mcmc
+        x0 = self.real_parameters * 1.1
+        mcmc = pints.RaoBlackwellACMC(x0)
+
+        # Configure
+        mcmc.set_target_acceptance_rate(0.3)
+        mcmc.set_initial_phase(True)
+
+        # Perform short run
+        rate = []
+        chain = []
+        for i in range(100):
+            x = mcmc.ask()
+            fx = self.log_posterior(x)
+            sample = mcmc.tell(fx)
+            if i == 20:
+                mcmc.set_initial_phase(False)
+            if i >= 50:
+                chain.append(sample)
+            rate.append(mcmc.acceptance_rate())
+            if np.all(sample == x):
+                self.assertEqual(mcmc.current_log_pdf(), fx)
+
+        chain = np.array(chain)
+        rate = np.array(rate)
+        self.assertEqual(chain.shape[0], 50)
+        self.assertEqual(chain.shape[1], len(x0))
+        self.assertEqual(rate.shape[0], 100)
+
+    def test_options(self):
+
+        # Test setting acceptance rate
+        x0 = self.real_parameters
+        mcmc = pints.RaoBlackwellACMC(x0)
+        self.assertNotEqual(mcmc.target_acceptance_rate(), 0.5)
+        mcmc.set_target_acceptance_rate(0.5)
+        self.assertEqual(mcmc.target_acceptance_rate(), 0.5)
+        mcmc.set_target_acceptance_rate(1)
+        self.assertRaises(ValueError, mcmc.set_target_acceptance_rate, 0)
+        self.assertRaises(ValueError, mcmc.set_target_acceptance_rate, -1e-6)
+        self.assertRaises(ValueError, mcmc.set_target_acceptance_rate, 1.00001)
+
+        # test hyperparameter setters and getters
+        self.assertEqual(mcmc.n_hyper_parameters(), 1)
+        self.assertRaises(ValueError, mcmc.set_hyper_parameters, [-0.1])
+        mcmc.set_hyper_parameters([0.3])
+        self.assertEqual(mcmc.eta(), 0.3)
+
+        self.assertEqual(mcmc.name(), 'Rao-Blackwell adaptive covariance MCMC')
+
+    def test_logging(self):
+
+        # Test logging includes name.
+        x = [self.real_parameters] * 3
+        mcmc = pints.MCMCController(
+            self.log_posterior, 3, x, method=pints.RaoBlackwellACMC)
+        mcmc.set_max_iterations(5)
+        with StreamCapture() as c:
+            mcmc.run()
+        text = c.text()
+        self.assertIn('Rao-Blackwell adaptive covariance MCMC', text)
+
+
+if __name__ == '__main__':
+    unittest.main()