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
elliptical_slice_sampler.py 文件源码
python
阅读 21
收藏 0
点赞 0
评论 0
评论列表
文章目录