Elliptical Slice Sampler

elliptical_slice_sampler.py

def ess_step(state, loglike_fn, prior_sampler, prev_loglike=None):
    """
    state: current state
    loglike_fn: loglikelihood function which accepts one argument: a
        parameter vector.
    prior_sampler: function (with no arguments) to sample from the prior,
        which must be a multivariate normal. i.e. the function returns a
        sample from a multivariate normal.
    prev_loglike (float): the result of `loglike_fn(state)`, which should be computed previously.
    """
    # 1. Choose ellipse.
    nu = prior_sampler()

    # 2. Compute log-likelihood threshold.
    u = np.random.rand()

    if prev_loglike is None:
        prev_loglike = loglike_fn(state)

    log_y = prev_loglike + np.log(u)

    # 3. Dran an initial proposal, also defining a bracket.
    two_pi = 2 * np.pi
    theta = np.random.uniform(0, two_pi)
    theta_min, theta_max = (theta - two_pi, theta)

    # 4. - 10.
    while True:
        cand = state * np.cos(theta) + nu * np.sin(theta)
        loglike_cand = loglike_fn(cand)
        if  loglike_cand > log_y:
            state = cand
            prev_loglike = loglike_cand
            return (state, prev_loglike)
        else:
            # Shrink the bracket and try a new point.
            if theta < 0:
                theta_min = theta
            else:
                theta_max = theta

            theta = np.random.uniform(theta_min, theta_max)

http://proceedings.mlr.press/v9/murray10a/murray10a.pdf