Source code for bumps.dream.metropolis

"""
MCMC step acceptance test.
"""

__all__ = ["metropolis", "metropolis_dr"]

import os

from numpy import array, cov, dot, errstate, exp, eye, minimum, sqrt, where
from numpy.linalg import cholesky, inv, norm

from . import util

BUMPS_TEMPERATURE = float(os.environ.get("BUMPS_TEMPERATURE", "1"))


def paccept(logp_old, logp_try):
    """
    Returns the probability of taking a metropolis step given two
    log density values.
    """
    return exp(minimum(logp_try - logp_old, 0) / BUMPS_TEMPERATURE)


[docs] def metropolis(xtry, logp_try, xold, logp_old, step_alpha): """ Metropolis rule for acceptance or rejection Generates the next generation, *newgen* from:: x_new[k] = x[k] if U > alpha = x_old[k] if U <= alpha where alpha is p/p_old and accept is U > alpha. Returns x_new, logp_new, alpha, accept """ with errstate(under="ignore"): alpha = paccept(logp_try=logp_try, logp_old=logp_old) alpha *= step_alpha accept = alpha > util.rng.rand(*alpha.shape) logp_new = where(accept, logp_try, logp_old) ## The following only works for vectors: # xnew = where(accept, xtry, xold) xnew = xtry + 0 for i, a in enumerate(accept): if not a: xnew[i] = xold[i] return xnew, logp_new, alpha, accept
def dr_step(x, scale): """ Delayed rejection step. """ # Compute the Cholesky Decomposition of X nchains, npars = x.shape r = (2.38 / sqrt(npars)) * cholesky(cov(x.T) + 1e-5 * eye(npars)) # Now do a delayed rejection step for each chain delta_x = dot(util.rng.randn(*x.shape), r) / scale # Generate ergodicity term eps = 1e-6 * util.rng.randn(*x.shape) # Update x_old with delta_x and eps; return x + delta_x + eps, r
[docs] def metropolis_dr(xtry, logp_try, x, logp, xold, logp_old, alpha12, R): """ Delayed rejection metropolis """ # Compute alpha32 (note we turned x and xtry around!) alpha32 = paccept(logp_try=logp, logp_old=logp_try) # Calculate alpha for each chain l2 = paccept(logp_try=logp_try, logp_old=logp_old) iR = inv(R) q1 = array( [exp(-0.5 * (norm(dot(x2 - x1, iR)) ** 2 - norm(dot(x1 - x0, iR)) ** 2)) for x0, x1, x2 in zip(xold, x, xtry)] ) alpha13 = l2 * q1 * (1 - alpha32) / (1 - alpha12) accept = alpha13 > util.rng.rand(*alpha13.shape) logp_new = where(accept, logp_try, logp) ## The following only works for vectors: # xnew = where(accept, xtry, x) xnew = xtry + 0 for i, a in enumerate(accept): if not a: xnew[i] = x[i] return xnew, logp_new, alpha13, accept