Skip to content

Commit b0a1280

Browse files
committed
Update computation of proposal logp ratio in SMC IMH kernel
1 parent 8f3636d commit b0a1280

File tree

1 file changed

+13
-11
lines changed

1 file changed

+13
-11
lines changed

pymc/smc/smc.py

Lines changed: 13 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -379,31 +379,33 @@ def tune(self):
379379
def mutate(self):
380380
"""Independent Metropolis-Hastings perturbation."""
381381
ac_ = np.empty((self.n_steps, self.draws))
382-
383-
cov = self.proposal_dist.cov
384382
log_R = np.log(self.rng.random((self.n_steps, self.draws)))
383+
384+
# The proposal is independent from the current point.
385+
# We have to take that into account to compute the Metropolis-Hastings acceptance
386+
# We first compute the logp of proposing a transition to the current points.
387+
# This variable is updated at the end of the loop with the entries from the accepted
388+
# transitions, which is equivalent to recomputing it in every iteration of the loop.
389+
backward_logp = self.proposal_dist.logpdf(self.tempered_posterior)
385390
for n_step in range(self.n_steps):
386-
# The proposal is independent from the current point.
387-
# We have to take that into account to compute the Metropolis-Hastings acceptance
388391
proposal = floatX(self.proposal_dist.rvs(size=self.draws, random_state=self.rng))
389392
proposal = proposal.reshape(len(proposal), -1)
390-
# To do that we compute the logp of moving to a new point
391-
forward = self.proposal_dist.logpdf(proposal)
392-
# And to going back from that new point
393-
backward = multivariate_normal(proposal.mean(axis=0), cov).logpdf(
394-
self.tempered_posterior
395-
)
393+
# We then compute the logp of proposing a transition to the new points
394+
forward_logp = self.proposal_dist.logpdf(proposal)
395+
396396
ll = np.array([self.likelihood_logp_func(prop) for prop in proposal])
397397
pl = np.array([self.prior_logp_func(prop) for prop in proposal])
398398
proposal_logp = pl + ll * self.beta
399399
accepted = log_R[n_step] < (
400-
(proposal_logp + backward) - (self.tempered_posterior_logp + forward)
400+
(proposal_logp + backward_logp) - (self.tempered_posterior_logp + forward_logp)
401401
)
402+
402403
ac_[n_step] = accepted
403404
self.tempered_posterior[accepted] = proposal[accepted]
404405
self.tempered_posterior_logp[accepted] = proposal_logp[accepted]
405406
self.prior_logp[accepted] = pl[accepted]
406407
self.likelihood_logp[accepted] = ll[accepted]
408+
backward_logp[accepted] = forward_logp[accepted]
407409

408410
self.acc_rate = np.mean(ac_)
409411

0 commit comments

Comments
 (0)