Reversible Jump MCMC examples

[Jeet Sukumaran]

Hi,

Are there any examples of setting up a RJ-MCMC in PyMC between models
of different dimensionality?

I've seen examples of model switching (with model index provided by a
Categorical variable), but the models have the same dimensionality.
I'm wondering if and how we might need to let the PyMC mcmc routines
know that the acceptance ratio needs to be modified when switching
dimensions.

Furthermore, assuming each model has at least some set of variables
which are not shared with other models, when currently "in" one model,
will PyMC's mcmc procedure needlessly update variables not used by
that model?
E.g. if M[0] has variables {x1, x2, x3}, and M[1] has variables {y1,
y2, y3}, and the model index variable C == 0, will the mcmc propose
changes to {y1, y2, y3}? Obviously, these will not change the log
probability of the data, given that this calculated with respect to
M[1]. This is not only inefficient, but will end up changing the state
of M[1] incorrectly and arbitrarily, right?

[Chris Fonnesbeck]

There is currently no built-in capability for RJMCMC in PyMC. In the sandbox directory, there is some RJMCMC code that was started a while back, but is incomplete. The biggest issue is calculating the Jacobian of the move across models every time a model jump occurs. We have not come up with a good way of generalizing this aspect of the algorithm.

[Jeet Sukumaran]

Thanks, Chris!

Yep, I can see it being challenging to come up with the general
approach to this. I guess a custom step sampler for the model indexing
stochastic variable would be the solution to this.

What about the issue of parent stochastics of the non-current model
being stepped through during the MCMC cycle? As they are not
contributing to the likelihood whether not the step is accepted is
arbitrary, and thus the sample will not be a real sample from the
posterior, right? I'm experimenting with a modification of the MCMC
object _loop method where it only steps stochastics that are the
parent of the current model. That is, instead of:

for step_method in self.step_methods:
if self.verbose > 2:
print_('Step method %s stepping' % step_method._id)
# Step the step method
step_method.step()

we have:

for step_method in self.step_method_dict[self.model_idx]:
step_method.step()
if self.model_idx == 0:
for step_method in
self.step_method_dict[self.geometric_p]:
step_method.step()
else:
for step_method in self.step_method_dict[self.poisson_mu]:
step_method.step()

(in this example, I am maintaining separate/independent stochastic
parameters for the Geometric and Poisson, instead of making one of
them a deterministic function of the other).


Is this a reasonable approach? Would I be adversely affecting some
other part of the MCMC machinery?