Blockwise Metropolis

Hi, I’m new to the pyMC community. I noticed that pm.sample(draws=10000, step=pm.Metropolis()) did a Metropolis within Gibbs sampling across all parameters. Is it possible to have pm.Metropolis() do blockwise sampling using a multivariate normal distribution? I’ve tried doing a compound step with disjoint parameter partitions using pm.Metropolis(). But this resulted in a Metropolis within Gibbs sampling across all parameters.

If you use a MvNormal distribution, the default sampler should attempt to update all elements of that random variable simultaneously.

That said, my experience is that using a multivariate Metropolis proposal without NUTS or HMC is usually pretty poor in performance and scales exponentially poorly with the dimension of the RV. The default sampler isn’t going to use blockwise conjugate updates for a multivariate normal, which is how I’ve often seen handwritten MCMC implementations do it.

Thank you for your response!

I’m using my own custom likelihood function. My code is as follows:

class LogLike(at.Op):

    itypes = [at.dvector]  # expects a vector of parameter values when called
    otypes = [at.dscalar]  # outputs a single scalar value (the log likelihood)

    def __init__(self, loglik, rtol, atol):
        Initialise the Op with various things that our log-likelihood function

            The log-likelihood

        # add inputs as class attributes
        self.likelihood = lambda sample: loglik(sample)

    def perform(self, node, inputs, outputs):
        # the method that is used when calling the Op
        (theta,) = inputs  # this will contain my variables
        # call the log-likelihood function
        logl = self.likelihood(theta)

logl = LogLike(loglik)
with pm.Model():
  variables = []
  for param_name in variable_name:
     variables.append(pm.Normal(param_name, mu=0., sigma=1.))
  theta = at.as_tensor_variable(variables)
  # use a Potential to "call" the Op and include it in the logp computation
  pm.Potential("likelihood", logl(theta))
  mh = pm.sample(draws=int(nsamples), step=pm.Metropolis())

To clarify, are you recommending I change

variables = []
for param_name in param_names:
   variables.append(pm.Normal(param_name, mu=0., sigma=1.))


variables = MvNormal("variables", mu=np.zeros(len(param_names)),cov(np.diag(ones(len(param_names))))?

Yes, that’s precisely right. Even though it’s the same underlying probability model, this way will force the sampler logic to treat the normal variates as a block.