elliptical_slice_sampler.py 文件源码

def elliptical_slice(xx, log_like_fn, prior_chol, prior_mean, *log_like_args, **sampler_args):
    cur_log_like = sampler_args.get('cur_log_like', None)
    angle_range = sampler_args.get('angle_range', 0)

    if cur_log_like is None:
        cur_log_like = log_like_fn(xx, *log_like_args)

    if np.isneginf(cur_log_like):
        raise Exception("Elliptical Slice Sampler: initial logprob is -inf for inputs %s" % xx)
    if np.isnan(cur_log_like):
        raise Exception("Elliptical Slice Sampler: initial logprob is NaN for inputs %s" % xx)

    nu = np.dot(prior_chol, npr.randn(xx.shape[0])) # don't bother adding mean here, would just subtract it at update step
    hh = np.log(npr.rand()) + cur_log_like  
    # log likelihood threshold -- LESS THAN THE INITIAL LOG LIKELIHOOD

    # Set up a bracket of angles and pick a first proposal.
    # "phi = (theta'-theta)" is a change in angle.
    if angle_range <= 0:
        # Bracket whole ellipse with both edges at first proposed point
        phi = npr.rand()*2*math.pi
        phi_min = phi - 2*math.pi
        phi_max = phi
    else:
        # Randomly center bracket on current point
        phi_min = -angle_range*npr.rand();
        phi_max = phi_min + angle_range;
        phi = npr.rand()*(phi_max - phi_min) + phi_min;

    # Slice sampling loop
    while True:
        # Compute xx for proposed angle difference 
        # and check if it's on the slice
        xx_prop = (xx-prior_mean)*np.cos(phi) + nu*np.sin(phi) + prior_mean

        cur_log_like = log_like_fn(xx_prop, *log_like_args)

        if cur_log_like > hh:
            # New point is on slice, ** EXIT LOOP **
            return xx_prop, cur_log_like

        # Shrink slice to rejected point
        if phi > 0:
            phi_max = phi
        elif phi < 0:
            phi_min = phi
        else:
            sys.stderr.write('Initial x: %s\n' % xx)
            # sys.stderr.write('initial log like = %f\n' % initial_log_like)
            sys.stderr.write('Proposed x: %s\n' % xx_prop)
            sys.stderr.write('ESS log lik = %f\n' % cur_log_like)
            raise Exception('BUG DETECTED: Shrunk to current position '
                            'and still not acceptable.');

        # Propose new angle difference
        phi = npr.rand()*(phi_max - phi_min) + phi_min