From 21a51fab266f66653f566132d9029f6eab11288f Mon Sep 17 00:00:00 2001 From: Gregor Boehl Date: Wed, 5 Oct 2022 15:09:38 +0200 Subject: [PATCH 01/14] add dime move --- src/emcee/moves/__init__.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/emcee/moves/__init__.py b/src/emcee/moves/__init__.py index 94670341..d99a6c03 100644 --- a/src/emcee/moves/__init__.py +++ b/src/emcee/moves/__init__.py @@ -2,6 +2,7 @@ from .de import DEMove from .de_snooker import DESnookerMove +from .dime import DIMEMove from .gaussian import GaussianMove from .kde import KDEMove from .mh import MHMove @@ -20,4 +21,5 @@ "KDEMove", "DEMove", "DESnookerMove", + "DIMEMove", ] From 4e1355926c29176b17760e68594069d61bc08803 Mon Sep 17 00:00:00 2001 From: Gregor Boehl Date: Mon, 10 Oct 2022 10:46:31 +0200 Subject: [PATCH 02/14] add dime move --- src/emcee/moves/dime.py | 140 +++++++++++++++++++++++ src/emcee/tests/integration/test_dime.py | 14 +++ 2 files changed, 154 insertions(+) create mode 100644 src/emcee/moves/dime.py create mode 100644 src/emcee/tests/integration/test_dime.py diff --git a/src/emcee/moves/dime.py b/src/emcee/moves/dime.py new file mode 100644 index 00000000..3f364760 --- /dev/null +++ b/src/emcee/moves/dime.py @@ -0,0 +1,140 @@ +# -*- coding: utf-8 -*- + +import numpy as np +from scipy.special import logsumexp +from scipy.stats import multivariate_t + +from .red_blue import RedBlueMove + +__all__ = ["DIMEMove"] + + +def mvt_sample(df, mean, cov, size, random): + """Sample from multivariate t distribution + + For reasons beyond my understanding, the results from random.multivariate_normal with non-identity covariance matrix are not reproducibel across architecture. Since scipy.stats.multivariate_t is based on numpy's multivariate_normal, the workaround is to crochet this manually. Advantage is that the scipy dependency drops out. + """ + + dim = len(mean) + + # draw samples + snorm = random.randn(size, dim) + chi2 = random.chisquare(df, size) / df + + # calculate sqrt of covariance + svd_cov = np.linalg.svd(cov * (df - 2) / df) + sqrt_cov = svd_cov[0] * np.sqrt(svd_cov[1]) @ svd_cov[2] + + return mean + snorm @ sqrt_cov / np.sqrt(chi2)[:, None] + + +class DIMEMove(RedBlueMove): + r"""A proposal using adaptive differential-independence mixture ensemble MCMC. + + This is the `Differential-Independence Mixture Ensemble proposal` as developed in `Ensemble MCMC Sampling for DSGE Models `_ (previousy ADEMC). + + Args: + sigma (Optional[float]): standard deviation of the Gaussian used to stretch the proposal vector. + gamma (Optional[float]): mean stretch factor for the proposal vector. By default, it is :math:`2.38 / \sqrt{2\,\mathrm{ndim}}` as recommended by `ter Braak (2006) `_. + aimh_prob (Optional[float]): probability to draw an adaptive independence Metropolis Hastings (AIMH) proposal. By default this is set to :math:`0.1`. + df_proposal_dist (Optional[float]): degrees of freedom of the multivariate t distribution used for AIMH proposals. Defaults to :math:`10`. + """ + + def __init__( + self, sigma=1.0e-5, gamma=None, aimh_prob=0.1, df_proposal_dist=10, **kwargs + ): + + self.sigma = sigma + self.g0 = gamma + self.aimh_prob = aimh_prob + self.dft = df_proposal_dist + + kwargs["nsplits"] = 1 + super(DIMEMove, self).__init__(**kwargs) + + def setup(self, coords): + # set some sane defaults + + nchain, npar = coords.shape + + if self.g0 is None: + # pure MAGIC + self.g0 = 2.38 / np.sqrt(2 * npar) + + if not hasattr(self, "prop_cov"): + # even more MAGIC + self.prop_cov = np.eye(npar) + self.prop_mean = np.zeros(npar) + self.accepted = np.ones(nchain, dtype=bool) + self.cumlweight = -np.inf + + def propose(self, model, state): + # wrap original propose to get the some info on the current state + self.lprobs = state.log_prob + state, accepted = super(DIMEMove, self).propose(model, state) + self.accepted = accepted + return state, accepted + + def get_proposal(self, x, dummy, random): + # actual proposal function + + nchain, npar = x.shape + + # differential evolution: draw the indices of the complementary chains + i0 = np.arange(nchain) + random.randint(1, nchain, size=nchain) + i1 = np.arange(nchain) + random.randint(1, nchain - 1, size=nchain) + i1[i1 >= i0] += 1 + # add small noise and calculate proposal + f = self.sigma * random.randn(nchain) + q = x + self.g0 * (x[i0 % nchain] - x[i1 % nchain]) + f[:, np.newaxis] + factors = np.zeros(nchain, dtype=np.float64) + + # log weight of current ensemble + lweight = logsumexp(self.lprobs) + \ + np.log(sum(self.accepted)) - np.log(nchain) + + # calculate stats for current ensemble + ncov = np.cov(x.T, ddof=1) + nmean = np.mean(x, axis=0) + + # update AIMH proposal distribution + newcumlweight = np.logaddexp(self.cumlweight, lweight) + self.prop_cov = ( + np.exp(self.cumlweight - newcumlweight) * self.prop_cov + + np.exp(lweight - newcumlweight) * ncov + ) + self.prop_mean = ( + np.exp(self.cumlweight - newcumlweight) * self.prop_mean + + np.exp(lweight - newcumlweight) * nmean + ) + self.cumlweight = newcumlweight + + # draw chains for AIMH sampling + xchnge = random.rand(nchain) <= self.aimh_prob + + # draw alternative candidates and calculate their proposal density + xcand = mvt_sample( + df=self.dft, + mean=self.prop_mean, + cov=self.prop_cov, + size=sum(xchnge), + random=random, + ) + lprop_old = multivariate_t.logpdf( + x[xchnge], + self.prop_mean, + self.prop_cov * (self.dft - 2) / self.dft, + df=self.dft, + ) + lprop_new = multivariate_t.logpdf( + xcand, + self.prop_mean, + self.prop_cov * (self.dft - 2) / self.dft, + df=self.dft, + ) + + # update proposals and factors + q[xchnge, :] = np.reshape(xcand, (-1, npar)) + factors[xchnge] = lprop_old - lprop_new + + return q, factors diff --git a/src/emcee/tests/integration/test_dime.py b/src/emcee/tests/integration/test_dime.py new file mode 100644 index 00000000..81ef3ca2 --- /dev/null +++ b/src/emcee/tests/integration/test_dime.py @@ -0,0 +1,14 @@ +# -*- coding: utf-8 -*- + +from emcee import moves + +from .test_proposal import _test_normal, _test_uniform + +__all__ = ["test_normal_dime", "test_uniform_de"] + + +def test_normal_dime(**kwargs): + _test_normal(moves.DIMEMove(), **kwargs) + +def test_uniform_dime(**kwargs): + _test_uniform(moves.DIMEMove(), **kwargs) From c1e49a3181b6d2c2b2f3f3e60a3f24085b68cc9d Mon Sep 17 00:00:00 2001 From: Gregor Boehl Date: Mon, 10 Oct 2022 10:47:44 +0200 Subject: [PATCH 03/14] update docstring --- src/emcee/moves/dime.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/emcee/moves/dime.py b/src/emcee/moves/dime.py index 3f364760..7b16219f 100644 --- a/src/emcee/moves/dime.py +++ b/src/emcee/moves/dime.py @@ -31,7 +31,7 @@ def mvt_sample(df, mean, cov, size, random): class DIMEMove(RedBlueMove): r"""A proposal using adaptive differential-independence mixture ensemble MCMC. - This is the `Differential-Independence Mixture Ensemble proposal` as developed in `Ensemble MCMC Sampling for DSGE Models `_ (previousy ADEMC). + This is the `Differential-Independence Mixture Ensemble proposal` as developed in `Ensemble MCMC Sampling for DSGE Models `_. Args: sigma (Optional[float]): standard deviation of the Gaussian used to stretch the proposal vector. From 023251226aea05e6879d31bd2d4f8eda8c383ca1 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 10 Oct 2022 09:22:45 +0000 Subject: [PATCH 04/14] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/emcee/moves/dime.py | 14 +++++++++++--- src/emcee/tests/integration/test_dime.py | 1 + 2 files changed, 12 insertions(+), 3 deletions(-) diff --git a/src/emcee/moves/dime.py b/src/emcee/moves/dime.py index 7b16219f..5e69cab0 100644 --- a/src/emcee/moves/dime.py +++ b/src/emcee/moves/dime.py @@ -41,7 +41,12 @@ class DIMEMove(RedBlueMove): """ def __init__( - self, sigma=1.0e-5, gamma=None, aimh_prob=0.1, df_proposal_dist=10, **kwargs + self, + sigma=1.0e-5, + gamma=None, + aimh_prob=0.1, + df_proposal_dist=10, + **kwargs ): self.sigma = sigma @@ -90,8 +95,11 @@ def get_proposal(self, x, dummy, random): factors = np.zeros(nchain, dtype=np.float64) # log weight of current ensemble - lweight = logsumexp(self.lprobs) + \ - np.log(sum(self.accepted)) - np.log(nchain) + lweight = ( + logsumexp(self.lprobs) + + np.log(sum(self.accepted)) + - np.log(nchain) + ) # calculate stats for current ensemble ncov = np.cov(x.T, ddof=1) diff --git a/src/emcee/tests/integration/test_dime.py b/src/emcee/tests/integration/test_dime.py index 81ef3ca2..49a3c2ff 100644 --- a/src/emcee/tests/integration/test_dime.py +++ b/src/emcee/tests/integration/test_dime.py @@ -10,5 +10,6 @@ def test_normal_dime(**kwargs): _test_normal(moves.DIMEMove(), **kwargs) + def test_uniform_dime(**kwargs): _test_uniform(moves.DIMEMove(), **kwargs) From b9e0b70b0b5161a6cd31ac81181c1c38f95b076e Mon Sep 17 00:00:00 2001 From: Gregor Boehl Date: Mon, 10 Oct 2022 11:35:35 +0200 Subject: [PATCH 05/14] skip test if no scipy --- src/emcee/tests/integration/test_dime.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/emcee/tests/integration/test_dime.py b/src/emcee/tests/integration/test_dime.py index 81ef3ca2..8351883b 100644 --- a/src/emcee/tests/integration/test_dime.py +++ b/src/emcee/tests/integration/test_dime.py @@ -1,5 +1,11 @@ # -*- coding: utf-8 -*- +try: + import scipy +except ImportError: + scipy = None +import pytest + from emcee import moves from .test_proposal import _test_normal, _test_uniform @@ -7,8 +13,10 @@ __all__ = ["test_normal_dime", "test_uniform_de"] +@pytest.mark.skipif(scipy is None, reason="scipy is not available") def test_normal_dime(**kwargs): _test_normal(moves.DIMEMove(), **kwargs) +@pytest.mark.skipif(scipy is None, reason="scipy is not available") def test_uniform_dime(**kwargs): _test_uniform(moves.DIMEMove(), **kwargs) From f895783940e5bc8e4bc6b010c5089400e92d64d5 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 10 Oct 2022 09:37:14 +0000 Subject: [PATCH 06/14] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/emcee/tests/integration/test_dime.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/emcee/tests/integration/test_dime.py b/src/emcee/tests/integration/test_dime.py index 8351883b..383efe71 100644 --- a/src/emcee/tests/integration/test_dime.py +++ b/src/emcee/tests/integration/test_dime.py @@ -17,6 +17,7 @@ def test_normal_dime(**kwargs): _test_normal(moves.DIMEMove(), **kwargs) + @pytest.mark.skipif(scipy is None, reason="scipy is not available") def test_uniform_dime(**kwargs): _test_uniform(moves.DIMEMove(), **kwargs) From 002e2d0c1c81f744b58182d6171554073f9415c8 Mon Sep 17 00:00:00 2001 From: Gregor Boehl Date: Mon, 10 Oct 2022 11:45:05 +0200 Subject: [PATCH 07/14] add import error if no scipy --- src/emcee/moves/dime.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/emcee/moves/dime.py b/src/emcee/moves/dime.py index 5e69cab0..92596b37 100644 --- a/src/emcee/moves/dime.py +++ b/src/emcee/moves/dime.py @@ -1,8 +1,11 @@ # -*- coding: utf-8 -*- import numpy as np -from scipy.special import logsumexp -from scipy.stats import multivariate_t +try: + from scipy.special import logsumexp + from scipy.stats import multivariate_t +except ImportError: + multivariate_t = None from .red_blue import RedBlueMove @@ -48,6 +51,10 @@ def __init__( df_proposal_dist=10, **kwargs ): + if multivariate_t is None: + raise ImportError( + "you need scipy.stats.multivariate_t and scipy.special.logsumexp to use the DIMEMove" + ) self.sigma = sigma self.g0 = gamma From 1af06cbf3d14fd6f32d9fc0a4ad8d698a01bd329 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 10 Oct 2022 09:45:17 +0000 Subject: [PATCH 08/14] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/emcee/moves/dime.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/emcee/moves/dime.py b/src/emcee/moves/dime.py index 92596b37..f061a132 100644 --- a/src/emcee/moves/dime.py +++ b/src/emcee/moves/dime.py @@ -1,6 +1,7 @@ # -*- coding: utf-8 -*- import numpy as np + try: from scipy.special import logsumexp from scipy.stats import multivariate_t From 52bb33a0ecb823d4e61f19736599e8b566055924 Mon Sep 17 00:00:00 2001 From: Gregor Boehl Date: Fri, 14 Oct 2022 02:56:25 +0200 Subject: [PATCH 09/14] some cleanup --- src/emcee/moves/dime.py | 55 ++++++++++++++++++++++------------------- 1 file changed, 29 insertions(+), 26 deletions(-) diff --git a/src/emcee/moves/dime.py b/src/emcee/moves/dime.py index f061a132..50d3f60a 100644 --- a/src/emcee/moves/dime.py +++ b/src/emcee/moves/dime.py @@ -82,32 +82,23 @@ def setup(self, coords): self.cumlweight = -np.inf def propose(self, model, state): - # wrap original propose to get the some info on the current state + """Wrap original propose to get the some info on the current state + """ + self.lprobs = state.log_prob state, accepted = super(DIMEMove, self).propose(model, state) self.accepted = accepted return state, accepted - def get_proposal(self, x, dummy, random): - # actual proposal function + def update_proposal_dist(self, x): + """Update proposal distribution with ensemble `x` + """ nchain, npar = x.shape - # differential evolution: draw the indices of the complementary chains - i0 = np.arange(nchain) + random.randint(1, nchain, size=nchain) - i1 = np.arange(nchain) + random.randint(1, nchain - 1, size=nchain) - i1[i1 >= i0] += 1 - # add small noise and calculate proposal - f = self.sigma * random.randn(nchain) - q = x + self.g0 * (x[i0 % nchain] - x[i1 % nchain]) + f[:, np.newaxis] - factors = np.zeros(nchain, dtype=np.float64) - # log weight of current ensemble - lweight = ( - logsumexp(self.lprobs) - + np.log(sum(self.accepted)) - - np.log(nchain) - ) + lweight = logsumexp(self.lprobs) + \ + np.log(sum(self.accepted)) - np.log(nchain) # calculate stats for current ensemble ncov = np.cov(x.T, ddof=1) @@ -125,6 +116,24 @@ def get_proposal(self, x, dummy, random): ) self.cumlweight = newcumlweight + def get_proposal(self, x, dummy, random): + """Actual proposal function + """ + + nchain, npar = x.shape + + # update AIMH proposal distribution + self.update_proposal_dist(x) + + # differential evolution: draw the indices of the complementary chains + i0 = np.arange(nchain) + random.randint(1, nchain, size=nchain) + i1 = np.arange(nchain) + random.randint(1, nchain - 1, size=nchain) + i1 += i1 >= i0 + # add small noise and calculate proposal + f = self.sigma * random.randn(nchain) + q = x + self.g0 * (x[i0 % nchain] - x[i1 % nchain]) + f[:, np.newaxis] + factors = np.zeros(nchain, dtype=np.float64) + # draw chains for AIMH sampling xchnge = random.rand(nchain) <= self.aimh_prob @@ -136,14 +145,8 @@ def get_proposal(self, x, dummy, random): size=sum(xchnge), random=random, ) - lprop_old = multivariate_t.logpdf( - x[xchnge], - self.prop_mean, - self.prop_cov * (self.dft - 2) / self.dft, - df=self.dft, - ) - lprop_new = multivariate_t.logpdf( - xcand, + lpropd = multivariate_t.logpdf( + np.vstack((x[None,xchnge],xcand[None])), self.prop_mean, self.prop_cov * (self.dft - 2) / self.dft, df=self.dft, @@ -151,6 +154,6 @@ def get_proposal(self, x, dummy, random): # update proposals and factors q[xchnge, :] = np.reshape(xcand, (-1, npar)) - factors[xchnge] = lprop_old - lprop_new + factors[xchnge] = lpropd[0] - lpropd[1] return q, factors From 3462623110ce609f45f9ac583afeec000916cebd Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Fri, 14 Oct 2022 00:57:33 +0000 Subject: [PATCH 10/14] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/emcee/moves/dime.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/emcee/moves/dime.py b/src/emcee/moves/dime.py index 50d3f60a..ed02aac3 100644 --- a/src/emcee/moves/dime.py +++ b/src/emcee/moves/dime.py @@ -82,8 +82,7 @@ def setup(self, coords): self.cumlweight = -np.inf def propose(self, model, state): - """Wrap original propose to get the some info on the current state - """ + """Wrap original propose to get the some info on the current state""" self.lprobs = state.log_prob state, accepted = super(DIMEMove, self).propose(model, state) @@ -91,14 +90,16 @@ def propose(self, model, state): return state, accepted def update_proposal_dist(self, x): - """Update proposal distribution with ensemble `x` - """ + """Update proposal distribution with ensemble `x`""" nchain, npar = x.shape # log weight of current ensemble - lweight = logsumexp(self.lprobs) + \ - np.log(sum(self.accepted)) - np.log(nchain) + lweight = ( + logsumexp(self.lprobs) + + np.log(sum(self.accepted)) + - np.log(nchain) + ) # calculate stats for current ensemble ncov = np.cov(x.T, ddof=1) @@ -117,8 +118,7 @@ def update_proposal_dist(self, x): self.cumlweight = newcumlweight def get_proposal(self, x, dummy, random): - """Actual proposal function - """ + """Actual proposal function""" nchain, npar = x.shape @@ -146,7 +146,7 @@ def get_proposal(self, x, dummy, random): random=random, ) lpropd = multivariate_t.logpdf( - np.vstack((x[None,xchnge],xcand[None])), + np.vstack((x[None, xchnge], xcand[None])), self.prop_mean, self.prop_cov * (self.dft - 2) / self.dft, df=self.dft, From 945f8ff9d9c8b8f5259697d173503b10d313995c Mon Sep 17 00:00:00 2001 From: Gregor Boehl Date: Fri, 14 Oct 2022 16:23:28 +0200 Subject: [PATCH 11/14] simplify --- src/emcee/moves/dime.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/emcee/moves/dime.py b/src/emcee/moves/dime.py index 50d3f60a..6b5ecae0 100644 --- a/src/emcee/moves/dime.py +++ b/src/emcee/moves/dime.py @@ -145,7 +145,7 @@ def get_proposal(self, x, dummy, random): size=sum(xchnge), random=random, ) - lpropd = multivariate_t.logpdf( + lpropdold, lpropdnew = multivariate_t.logpdf( np.vstack((x[None,xchnge],xcand[None])), self.prop_mean, self.prop_cov * (self.dft - 2) / self.dft, @@ -154,6 +154,6 @@ def get_proposal(self, x, dummy, random): # update proposals and factors q[xchnge, :] = np.reshape(xcand, (-1, npar)) - factors[xchnge] = lpropd[0] - lpropd[1] + factors[xchnge] = lpropdold - lpropdnew return q, factors From 9c811d7b73047460fc2528794f6c77e83d6947f7 Mon Sep 17 00:00:00 2001 From: Gregor Boehl Date: Fri, 14 Oct 2022 16:25:46 +0200 Subject: [PATCH 12/14] simplify --- src/emcee/moves/dime.py | 1 + 1 file changed, 1 insertion(+) diff --git a/src/emcee/moves/dime.py b/src/emcee/moves/dime.py index 56d47564..320f4ab6 100644 --- a/src/emcee/moves/dime.py +++ b/src/emcee/moves/dime.py @@ -146,6 +146,7 @@ def get_proposal(self, x, dummy, random): random=random, ) lpropdold, lpropdnew = multivariate_t.logpdf( + np.vstack((x[None, xchnge], xcand[None])), self.prop_mean, self.prop_cov * (self.dft - 2) / self.dft, df=self.dft, From 6399a33582d0d021d988a2773bf1315d478c6580 Mon Sep 17 00:00:00 2001 From: Gregor Boehl Date: Sun, 2 Jun 2024 07:39:48 +0200 Subject: [PATCH 13/14] update & sync --- src/emcee/moves/dime.py | 64 +++++++++++++++++++++++------------------ 1 file changed, 36 insertions(+), 28 deletions(-) diff --git a/src/emcee/moves/dime.py b/src/emcee/moves/dime.py index 320f4ab6..3e25b80b 100644 --- a/src/emcee/moves/dime.py +++ b/src/emcee/moves/dime.py @@ -16,7 +16,7 @@ def mvt_sample(df, mean, cov, size, random): """Sample from multivariate t distribution - For reasons beyond my understanding, the results from random.multivariate_normal with non-identity covariance matrix are not reproducibel across architecture. Since scipy.stats.multivariate_t is based on numpy's multivariate_normal, the workaround is to crochet this manually. Advantage is that the scipy dependency drops out. + The results from random.multivariate_normal with non-identity covariance matrix are not reproducibel across OS architecture. Since scipy.stats.multivariate_t is based on numpy's multivariate_normal, the workaround is to provide this manually. """ dim = len(mean) @@ -37,20 +37,22 @@ class DIMEMove(RedBlueMove): This is the `Differential-Independence Mixture Ensemble proposal` as developed in `Ensemble MCMC Sampling for DSGE Models `_. - Args: - sigma (Optional[float]): standard deviation of the Gaussian used to stretch the proposal vector. - gamma (Optional[float]): mean stretch factor for the proposal vector. By default, it is :math:`2.38 / \sqrt{2\,\mathrm{ndim}}` as recommended by `ter Braak (2006) `_. - aimh_prob (Optional[float]): probability to draw an adaptive independence Metropolis Hastings (AIMH) proposal. By default this is set to :math:`0.1`. - df_proposal_dist (Optional[float]): degrees of freedom of the multivariate t distribution used for AIMH proposals. Defaults to :math:`10`. + Parameters + ---------- + sigma : float, optional + standard deviation of the Gaussian used to stretch the proposal vector. + gamma : float, optional + mean stretch factor for the proposal vector. By default, it is :math:`2.38 / \sqrt{2\,\mathrm{ndim}}` as recommended by `ter Braak (2006) `_. + aimh_prob : float, optional + probability to draw an adaptive independence Metropolis Hastings (AIMH) proposal. By default this is set to :math:`0.1`. + df_proposal_dist : float, optional + degrees of freedom of the multivariate t distribution used for AIMH proposals. Defaults to :math:`10`. + rho : float, optional + decay parameter for the aimh proposal mean and covariances. Defaults to :math:`0.999`. """ def __init__( - self, - sigma=1.0e-5, - gamma=None, - aimh_prob=0.1, - df_proposal_dist=10, - **kwargs + self, sigma=1.0e-5, gamma=None, aimh_prob=0.1, df_proposal_dist=10, rho=.999, **kwargs ): if multivariate_t is None: raise ImportError( @@ -59,10 +61,10 @@ def __init__( self.sigma = sigma self.g0 = gamma + self.decay = rho self.aimh_prob = aimh_prob self.dft = df_proposal_dist - kwargs["nsplits"] = 1 super(DIMEMove, self).__init__(**kwargs) def setup(self, coords): @@ -80,9 +82,13 @@ def setup(self, coords): self.prop_mean = np.zeros(npar) self.accepted = np.ones(nchain, dtype=bool) self.cumlweight = -np.inf + else: + # update AIMH proposal distribution + self.update_proposal_dist(coords) def propose(self, model, state): - """Wrap original propose to get the some info on the current state""" + """Wrap original propose to get the some info on the current state + """ self.lprobs = state.log_prob state, accepted = super(DIMEMove, self).propose(model, state) @@ -90,16 +96,17 @@ def propose(self, model, state): return state, accepted def update_proposal_dist(self, x): - """Update proposal distribution with ensemble `x`""" + """Update proposal distribution with ensemble `x` + """ nchain, npar = x.shape # log weight of current ensemble - lweight = ( - logsumexp(self.lprobs) - + np.log(sum(self.accepted)) - - np.log(nchain) - ) + if self.accepted.any(): + lweight = logsumexp(self.lprobs) + \ + np.log(sum(self.accepted)) - np.log(nchain) + else: + lweight = -np.inf # calculate stats for current ensemble ncov = np.cov(x.T, ddof=1) @@ -115,15 +122,16 @@ def update_proposal_dist(self, x): np.exp(self.cumlweight - newcumlweight) * self.prop_mean + np.exp(lweight - newcumlweight) * nmean ) + # self.cumlweight = newcumlweight + np.log(self.decay) self.cumlweight = newcumlweight - def get_proposal(self, x, dummy, random): - """Actual proposal function""" + def get_proposal(self, x, xref, random): + """Actual proposal function + """ + xref = xref[0] nchain, npar = x.shape - - # update AIMH proposal distribution - self.update_proposal_dist(x) + nref, _ = xref.shape # differential evolution: draw the indices of the complementary chains i0 = np.arange(nchain) + random.randint(1, nchain, size=nchain) @@ -131,7 +139,7 @@ def get_proposal(self, x, dummy, random): i1 += i1 >= i0 # add small noise and calculate proposal f = self.sigma * random.randn(nchain) - q = x + self.g0 * (x[i0 % nchain] - x[i1 % nchain]) + f[:, np.newaxis] + q = x + self.g0 * (xref[i0 % nref] - xref[i1 % nref]) + f[:, None] factors = np.zeros(nchain, dtype=np.float64) # draw chains for AIMH sampling @@ -145,7 +153,7 @@ def get_proposal(self, x, dummy, random): size=sum(xchnge), random=random, ) - lpropdold, lpropdnew = multivariate_t.logpdf( + lprop_old, lprop_new = multivariate_t.logpdf( np.vstack((x[None, xchnge], xcand[None])), self.prop_mean, self.prop_cov * (self.dft - 2) / self.dft, @@ -154,6 +162,6 @@ def get_proposal(self, x, dummy, random): # update proposals and factors q[xchnge, :] = np.reshape(xcand, (-1, npar)) - factors[xchnge] = lpropdold - lpropdnew + factors[xchnge] = lprop_old - lprop_new return q, factors From b47484944688d89d0db0c6082b2350dc2bb8efc5 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Sun, 2 Jun 2024 05:39:19 +0000 Subject: [PATCH 14/14] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/emcee/moves/dime.py | 24 +++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/src/emcee/moves/dime.py b/src/emcee/moves/dime.py index 3e25b80b..3abacec2 100644 --- a/src/emcee/moves/dime.py +++ b/src/emcee/moves/dime.py @@ -52,7 +52,13 @@ class DIMEMove(RedBlueMove): """ def __init__( - self, sigma=1.0e-5, gamma=None, aimh_prob=0.1, df_proposal_dist=10, rho=.999, **kwargs + self, + sigma=1.0e-5, + gamma=None, + aimh_prob=0.1, + df_proposal_dist=10, + rho=0.999, + **kwargs ): if multivariate_t is None: raise ImportError( @@ -87,8 +93,7 @@ def setup(self, coords): self.update_proposal_dist(coords) def propose(self, model, state): - """Wrap original propose to get the some info on the current state - """ + """Wrap original propose to get the some info on the current state""" self.lprobs = state.log_prob state, accepted = super(DIMEMove, self).propose(model, state) @@ -96,15 +101,17 @@ def propose(self, model, state): return state, accepted def update_proposal_dist(self, x): - """Update proposal distribution with ensemble `x` - """ + """Update proposal distribution with ensemble `x`""" nchain, npar = x.shape # log weight of current ensemble if self.accepted.any(): - lweight = logsumexp(self.lprobs) + \ - np.log(sum(self.accepted)) - np.log(nchain) + lweight = ( + logsumexp(self.lprobs) + + np.log(sum(self.accepted)) + - np.log(nchain) + ) else: lweight = -np.inf @@ -126,8 +133,7 @@ def update_proposal_dist(self, x): self.cumlweight = newcumlweight def get_proposal(self, x, xref, random): - """Actual proposal function - """ + """Actual proposal function""" xref = xref[0] nchain, npar = x.shape